From f61c0b5d03660a20d67ca75bee0b159db8e5dfa5 Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Fri, 27 Oct 2017 14:09:02 +0200 Subject: [PATCH] Very early version of MR solver --- lib/algorithms/Algorithms.h | 1 + lib/algorithms/iterative/MinimalResidual.h | 261 ++++++--------------- 2 files changed, 67 insertions(+), 195 deletions(-) diff --git a/lib/algorithms/Algorithms.h b/lib/algorithms/Algorithms.h index af83df67..503092db 100644 --- a/lib/algorithms/Algorithms.h +++ b/lib/algorithms/Algorithms.h @@ -47,6 +47,7 @@ Author: Peter Boyle #include #include #include +#include #include #include #include diff --git a/lib/algorithms/iterative/MinimalResidual.h b/lib/algorithms/iterative/MinimalResidual.h index a5104e03..3229b408 100644 --- a/lib/algorithms/iterative/MinimalResidual.h +++ b/lib/algorithms/iterative/MinimalResidual.h @@ -38,177 +38,24 @@ namespace Grid { // single input vec, single output vec. ///////////////////////////////////////////////////////////// -template -class MinimalResidual : public OperatorFunction { +template class MinimalResidual : public OperatorFunction { public: - bool ErrorOnNoConverge; // throw an assert when the MR fails to converge. - // Defaults true. - RealD Tolerance; + bool ErrorOnNoConverge; // throw an assert when the MR fails to converge. + // Defaults true. + RealD Tolerance; 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) - : Tolerance(tol), - MaxIterations(maxit), - ErrorOnNoConverge(err_on_no_conv){}; - - void operator()(LinearOperatorBase &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; - } + : Tolerance(tol), + MaxIterations(maxit), + ErrorOnNoConverge(err_on_no_conv){}; //! Minimal-residual (MR) algorithm for a generic Linear Operator /*! \ingroup invert * 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 * @@ -256,15 +103,13 @@ class MinimalResidual : public OperatorFunction { * @{ */ - // TODO: figure out what isign from chroma is supposed to do - void tmpImplFromChroma(LinearOperatorBase &Linop, const Field &src, - Field &psi) { + void operator()(LinearOperatorBase &Linop, const Field &src, Field &psi) { + psi.checkerboard = src.checkerboard; conformable(psi, src); Complex a, c; - Complex c; - RealD d; + RealD d; Field Mr(src); Field r(src); @@ -274,72 +119,98 @@ class MinimalResidual : public OperatorFunction { assert(std::isnan(guess) == 0); 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 := 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); - 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; } 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; - while( (k < MaxIterations) && (cp > rsd_sq) ) - { + while((k < MaxIterations) && (cp > rsd_sq)) { ++k; /* 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); 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); 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; - res.resid = sqrt(cp); - swatch.stop(); - std::cout << "InvMR: k = " << k << " cp = " << cp << endl; + // res.resid = sqrt(cp); + std::cout << "InvMR: k = " << k << " cp = " << cp << std::endl; // flopcount.report("invmr", swatch.getTimeInSeconds()); + std::cout << GridLogMessage << "MinimalResidual Converged on iteration " << k << std::endl; + std::cout << GridLogMessage << "\tComputed residual " << sqrt(cp / ssq)<