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

Very early version of MR solver

Browse Source
This commit is contained in:
Daniel Richtmann 2017-10-27 14:09:02 +02:00
parent 074db32e54
commit f61c0b5d03
No known key found for this signature in database
GPG Key ID: B33C490AF0772057
2 changed files with 67 additions and 195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
lib/algorithms/Algorithms.h
View File

@ -47,6 +47,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <Grid/algorithms/iterative/ConjugateGradientMixedPrec.h> #include <Grid/algorithms/iterative/ConjugateGradientMixedPrec.h>
#include <Grid/algorithms/iterative/BlockConjugateGradient.h> #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
#include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h> #include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h>
#include <Grid/algorithms/iterative/MinimalResidual.h>
#include <Grid/algorithms/iterative/GeneralisedMinimalResidual.h> #include <Grid/algorithms/iterative/GeneralisedMinimalResidual.h>
#include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h> #include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
#include <Grid/algorithms/CoarsenedMatrix.h> #include <Grid/algorithms/CoarsenedMatrix.h>

261
lib/algorithms/iterative/MinimalResidual.h
View File

@ -38,177 +38,24 @@ namespace Grid {
// single input vec, single output vec. // single input vec, single output vec.
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
template <class Field> template<class Field> class MinimalResidual : public OperatorFunction<Field> {
class MinimalResidual : public OperatorFunction<Field> {
public: public:
bool ErrorOnNoConverge; // throw an assert when the MR fails to converge. bool ErrorOnNoConverge; // throw an assert when the MR fails to converge.
// Defaults true. // Defaults true.
RealD Tolerance; RealD Tolerance;
Integer MaxIterations; Integer MaxIterations;
Integer IterationsToComplete; //Number of iterations the MR took to finish. Filled in upon completion Integer IterationsToComplete; // Number of iterations the MR took to finish. Filled in upon completion
MinimalResidual(RealD tol, Integer maxit, bool err_on_no_conv = true) MinimalResidual(RealD tol, Integer maxit, bool err_on_no_conv = true)
: Tolerance(tol), : Tolerance(tol),
MaxIterations(maxit), MaxIterations(maxit),
ErrorOnNoConverge(err_on_no_conv){}; ErrorOnNoConverge(err_on_no_conv){};
void operator()(LinearOperatorBase<Field> &Linop, const Field &src,
Field &psi) {
psi.checkerboard = src.checkerboard; // Check
conformable(psi, src);
/////
RealD cp, c, a, d, b, ssq, qq, b_pred;
Field p(src);
Field mmp(src);
Field r(src);
// Initial residual computation & set up
RealD guess = norm2(psi);
assert(std::isnan(guess) == 0);
/////
Field p {src};
Field matrixTimesPsi {src};
Field r {src};
RealD alpha {};
// Initial residual computation & set up
RealD guess = norm2(psi);
assert(std::isnan(guess) == 0);
Linop.HermOp(psi, matrixTimesPsi);
r = src - matrixTimesPsi;
Linop.HermOp(r, p);
alpha = innerProduct(p,r) / innerProduct(p,p);
psi = psi + alpha * r;
r = r - alpha * p;
Linop.HermOp(r, p);
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// RealD cp, c, a, d, b, ssq, qq, b_pred;
Field p(src);
Field matrixTimesPsi(src);
// Field r(src);
// Initial residual computation & set up
RealD guess = norm2(psi);
assert(std::isnan(guess) == 0);
Linop.HermOpAndNorm(psi, matrixTimesPsi, d, b);
r = src - matrixTimesPsi;
p = matrixTimesPsi;
a = norm2(p);
cp = a;
ssq = norm2(src);
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: guess " << guess << std::endl;
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: src " << ssq << std::endl;
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: mp " << d << std::endl;
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: matrixTimesPsi " << b << std::endl;
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: cp,r " << cp << std::endl;
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: p " << a << std::endl;
RealD rsq = Tolerance * Tolerance * ssq;
// Check if guess is really REALLY good :)
if (cp <= rsq) {
return;
}
std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: k=0 residual " << cp << " target " << rsq
<< std::endl;
GridStopWatch LinalgTimer;
GridStopWatch MatrixTimer;
GridStopWatch SolverTimer;
SolverTimer.Start();
int k;
for (k = 1; k <= MaxIterations; k++) {
c = cp;
MatrixTimer.Start();
Linop.HermOpAndNorm(p, matrixTimesPsi, d, qq);
MatrixTimer.Stop();
LinalgTimer.Start();
// RealD qqck = norm2(matrixTimesPsi);
// ComplexD dck = innerProduct(p,matrixTimesPsi);
a = c / d;
b_pred = a * (a * qq - d) / c;
cp = axpy_norm(r, -a, matrixTimesPsi, r);
b = cp / c;
// Fuse these loops ; should be really easy
psi = a * p + psi;
p = p * b + r;
LinalgTimer.Stop();
std::cout << GridLogIterative << "MinimalResidual: Iteration " << k
<< " residual " << cp << " target " << rsq << std::endl;
// Stopping condition
if (cp <= rsq) {
SolverTimer.Stop();
Linop.HermOpAndNorm(psi, matrixTimesPsi, d, qq);
p = matrixTimesPsi - src;
RealD matrixTimesPsiNorm = sqrt(norm2(matrixTimesPsi));
RealD psinorm = sqrt(norm2(psi));
RealD srcnorm = sqrt(norm2(src));
RealD resnorm = sqrt(norm2(p));
RealD true_residual = resnorm / srcnorm;
std::cout << GridLogMessage
<< "MinimalResidual: Converged on iteration " << k << std::endl;
std::cout << GridLogMessage << "Computed residual " << sqrt(cp / ssq)
<< " true residual " << true_residual << " target "
<< Tolerance << std::endl;
std::cout << GridLogMessage << "Time elapsed: Iterations "
<< SolverTimer.Elapsed() << " Matrix "
<< MatrixTimer.Elapsed() << " Linalg "
<< LinalgTimer.Elapsed();
std::cout << std::endl;
if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
IterationsToComplete = k;
return;
}
}
std::cout << GridLogMessage << "MinimalResidual did NOT converge"
<< std::endl;
if (ErrorOnNoConverge) assert(0);
IterationsToComplete = k;
}
//! Minimal-residual (MR) algorithm for a generic Linear Operator //! Minimal-residual (MR) algorithm for a generic Linear Operator
/*! \ingroup invert /*! \ingroup invert
* This subroutine uses the Minimal Residual (MR) algorithm to determine * This subroutine uses the Minimal Residual (MR) algorithm to determine
* the solution of the set of linear equations. Here we allow M to be nonhermitian. * the solution of the set of linear equations. Here we allow M to be
nonhermitian.
* *
* M . Psi = src * M . Psi = src
* *
@ -256,15 +103,13 @@ class MinimalResidual : public OperatorFunction<Field> {
* @{ * @{
*/ */
// TODO: figure out what isign from chroma is supposed to do void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
void tmpImplFromChroma(LinearOperatorBase<Field> &Linop, const Field &src,
Field &psi) {
psi.checkerboard = src.checkerboard; psi.checkerboard = src.checkerboard;
conformable(psi, src); conformable(psi, src);
Complex a, c; Complex a, c;
Complex c; RealD d;
RealD d;
Field Mr(src); Field Mr(src);
Field r(src); Field r(src);
@ -274,72 +119,98 @@ class MinimalResidual : public OperatorFunction<Field> {
assert(std::isnan(guess) == 0); assert(std::isnan(guess) == 0);
RealD ssq = norm2(src); // flopcount.addSiteFlops(4*Nc*Ns,s); // stands for "source squared" RealD ssq = norm2(src); // flopcount.addSiteFlops(4*Nc*Ns,s); // stands for "source squared"
RealD rsd_sq = Tolerance * Tolerance * ssq; // flopcount.addSiteFlops(4*Nc*Ns,s); // stands for "residual squared" RealD rsd_sq = Tolerance * Tolerance * ssq; // flopcount.addSiteFlops(4*Nc*Ns,s); //
// stands for "residual squared"
/* r[0] := src - M . Psi[0] */ /* r[0] := src - M . Psi[0] */
/* r := M . Psi */ /* r := M . Psi */
M(Mr, psi, isign); // flopcount.addFlops(M.nFlops()); // M(Mr, psi, isign); // flopcount.addFlops(M.nFlops());
Linop.Op(psi, Mr); // flopcount.addFlops(M.nFlops());
r = src - Mr; // flopcount.addSiteFlops(2*Nc*Ns,s); r = src - Mr; // flopcount.addSiteFlops(2*Nc*Ns,s);
RealD cp = norm2(r); /* Cp = |r[0]|^2 */ /* 2 Nc Ns flops */ // flopcount.addSiteFlops(4*Nc*Ns, s); RealD cp = norm2(r); /* Cp = |r[0]|^2 */
/* 2 Nc Ns flops */ // flopcount.addSiteFlops(4*Nc*Ns, s);
// auto cp = norm2(r); /* Cp = |r[0]|^2 */ /* 2 Nc Ns flops */ //
// flopcount.addSiteFlops(4*Nc*Ns, s);
if (cp <= rsd_sq) { /* IF |r[0]| <= Tolerance|src| THEN RETURN; */ if(cp <= rsd_sq) { /* IF |r[0]| <= Tolerance|src| THEN RETURN; */
return; return;
} }
std::cout << GridLogIterative << std::setprecision(4) std::cout << GridLogIterative << std::setprecision(4)
<< "MinimalResidual: k=0 residual " << cp << " target " << rsq_sq << std::endl; << "MinimalResidual: k=0 residual " << cp << " target " << rsd_sq << std::endl;
/* FOR k FROM 1 TO MaxIterations DO */ GridStopWatch LinalgTimer;
GridStopWatch MatrixTimer;
GridStopWatch SolverTimer;
SolverTimer.Start();
auto k = 0; auto k = 0;
while( (k < MaxIterations) && (cp > rsd_sq) ) while((k < MaxIterations) && (cp > rsd_sq)) {
{
++k; ++k;
/* a[k-1] := < M.r[k-1], r[k-1] >/ < M.r[k-1], M.r[k-1] > ; */ /* a[k-1] := < M.r[k-1], r[k-1] >/ < M.r[k-1], M.r[k-1] > ; */
M(Mr, r, isign); /* Mr = M * r */ // flopcount.addFlops(M.nFlops()); MatrixTimer.Start();
// M(Mr, r, isign); /* Mr = M * r */ // flopcount.addFlops(M.nFlops());
Linop.Op(r, Mr); /* Mr = M * r */ // flopcount.addFlops(M.nFlops());
MatrixTimer.Stop();
LinalgTimer.Start();
c = innerProduct(Mr, r); /* c = < M.r, r > */ // flopcount.addSiteFlops(4*Nc*Ns,s); c = innerProduct(Mr, r); /* c = < M.r, r > */ // flopcount.addSiteFlops(4*Nc*Ns,s);
d = norm2(Mr); /* d = | M.r | ** 2 */ // flopcount.addSiteFlops(4*Nc*Ns,s); d = norm2(Mr); /* d = | M.r | ** 2 */ // flopcount.addSiteFlops(4*Nc*Ns,s);
a = c / d; /* a = c / d */ a = c / d;
a = a * MRovpar; /* a[k-1] *= MRovpar ; */ // a = a * MRovpar; /* a[k-1] *= MRovpar ; */
psi = psi + r * a; /* Psi[k] += a[k-1] r[k-1] ; */ // flopcount.addSiteFlops(4*Nc*Ns,s);
psi = psi + r * a; /* Psi[k] += a[k-1] r[k-1] ; */ // flopcount.addSiteFlops(4*Nc*Ns,s);
r = r - Mr * a; /* r[k] -= a[k-1] M . r[k-1] ; */ // flopcount.addSiteFlops(4*Nc*Ns,s); r = r - Mr * a; /* r[k] -= a[k-1] M . r[k-1] ; */ // flopcount.addSiteFlops(4*Nc*Ns,s);
cp = norm2(r); /* cp = | r[k] |**2 */ // flopcount.addSiteFlops(4*Nc*Ns,s); cp = norm2(r); /* cp = | r[k] |**2 */ // flopcount.addSiteFlops(4*Nc*Ns,s);
// std::cout << "InvMR: k = " << k << " cp = " << cp << endl; LinalgTimer.Stop();
std::cout << GridLogIterative << "MinimalResidual: Iteration " << k
<< " residual " << cp << " target " << rsd_sq << std::endl;
} }
SolverTimer.Stop();
IterationsToComplete = k; IterationsToComplete = k;
res.resid = sqrt(cp); // res.resid = sqrt(cp);
swatch.stop(); std::cout << "InvMR: k = " << k << " cp = " << cp << std::endl;
std::cout << "InvMR: k = " << k << " cp = " << cp << endl;
// flopcount.report("invmr", swatch.getTimeInSeconds()); // flopcount.report("invmr", swatch.getTimeInSeconds());
std::cout << GridLogMessage << "MinimalResidual Converged on iteration " << k << std::endl;
std::cout << GridLogMessage << "\tComputed residual " << sqrt(cp / ssq)<<std::endl;
// std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
// std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
std::cout << GridLogMessage << "Time breakdown "<<std::endl;
std::cout << GridLogMessage << "\tElapsed " << SolverTimer.Elapsed() <<std::endl;
std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
std::cout << GridLogMessage << "\tLinalg " << LinalgTimer.Elapsed() <<std::endl;
// Compute the actual residual // Compute the actual residual
{ {
M(Mr, psi, isign); // M(Mr, psi, isign);
RealD actual_res = norm2(src- Mr); Linop.Op(psi, Mr);
res.resid = sqrt(actual_res); Field tmp = src - Mr;
// RealD actual_res = norm2(src-Mr);
RealD actual_res = norm2(tmp);
// res.resid = sqrt(actual_res);
} }
if ( IterationsToComplete == MaxIterations ) if(IterationsToComplete == MaxIterations)
std::cerr << "Nonconvergence Warning" << endl; std::cerr << "Nonconvergence Warning" << std::endl;
END_CODE();
return res;
// return res;
} }
}; };
} } // namespace Grid
#endif #endif