diff --git a/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h b/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h index 94193484..ce82ed4d 100644 --- a/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h @@ -29,30 +29,13 @@ directory #ifndef GRID_FLEXIBLE_GENERALISED_MINIMAL_RESIDUAL_H #define GRID_FLEXIBLE_GENERALISED_MINIMAL_RESIDUAL_H -// from Y. Saad - Iterative Methods for Sparse Linear Systems, PP 172 -// Compute r0 = b − Ax0 , β := ||r0||2 , and v1 := r0 /β -// For j = 1, 2, ..., m Do: -// Compute wj := Avj -// For i = 1, ..., j Do: -// hij := (wj , vi) -// wj := wj − hij vi -// EndDo -// hj+1,j = ||wj||2 . If hj+1,j = 0 set m := j and go to HERE -// vj+1 = wj /hj+1,j -// EndDo -// Define the (m + 1) × m Hessenberg matrix H̄m = {hij}1≤i≤m+1,1≤j≤m. [HERE] -// Compute ym the minimizer of ||βe1 − H̄m y||2 and xm = x0 + Vm ym. -/////////////////////////////////////////////////////////////////////////////////////////////////////// - -// want to solve Ax = b -> A = LinOp, psi = x, b = src - namespace Grid { template class FlexibleGeneralisedMinimalResidual : public OperatorFunction { public: bool ErrorOnNoConverge; // Throw an assert when FGMRES fails to converge, - // defaults to True. + // defaults to true RealD Tolerance; @@ -264,8 +247,6 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction { y[i] = y[i] / H(i, i); } - // TODO: Use axpys or similar for these - // TODO: Fix the condition if (true) { for (int i = 0; i <= iter; i++) psi = psi + v[i] * y[i]; @@ -280,5 +261,3 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction { }; } #endif - -// I took the version with p->kind == left from the WMG code base here. TODO: recheck if we need left or right diff --git a/lib/algorithms/iterative/GeneralisedMinimalResidual.h b/lib/algorithms/iterative/GeneralisedMinimalResidual.h index 27ab6caf..0f7bf155 100644 --- a/lib/algorithms/iterative/GeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/GeneralisedMinimalResidual.h @@ -35,7 +35,7 @@ template class GeneralisedMinimalResidual : public OperatorFunction { public: bool ErrorOnNoConverge; // Throw an assert when GMRES fails to converge, - // defaults to True. + // defaults to true RealD Tolerance; @@ -231,8 +231,6 @@ class GeneralisedMinimalResidual : public OperatorFunction { y[i] = y[i] / H(i, i); } - // TODO: Use axpys or similar for these - // TODO: Fix the condition if (true) { for (int i = 0; i <= iter; i++) psi = psi + v[i] * y[i];