From c6cbe533ea9690aace4805a2cce6e8228ef5feb3 Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Thu, 9 Nov 2017 17:12:54 +0100 Subject: [PATCH] Set everything up for the implementation of CAGMRES The current implementation is the exact same code as normal GMRES. This commit only sets up the "framework" for the implementation of CAGMRES, i.e., a test and an include in the algorithms header file. --- lib/algorithms/Algorithms.h | 1 + ...cationAvoidingGeneralisedMinimalResidual.h | 247 ++++++++++++++++++ tests/solver/Test_wilson_cagmres_unprec.cc | 65 +++++ 3 files changed, 313 insertions(+) create mode 100644 lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h create mode 100644 tests/solver/Test_wilson_cagmres_unprec.cc diff --git a/lib/algorithms/Algorithms.h b/lib/algorithms/Algorithms.h index 3368dce8..1517c0be 100644 --- a/lib/algorithms/Algorithms.h +++ b/lib/algorithms/Algorithms.h @@ -49,6 +49,7 @@ Author: Peter Boyle #include #include #include +#include #include #include #include diff --git a/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h b/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h new file mode 100644 index 00000000..2e574639 --- /dev/null +++ b/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h @@ -0,0 +1,247 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h + +Copyright (C) 2015 + +Author: Daniel Richtmann + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ +#ifndef GRID_COMMUNICATION_AVOIDING_GENERALISED_MINIMAL_RESIDUAL_H +#define GRID_COMMUNICATION_AVOIDING_GENERALISED_MINIMAL_RESIDUAL_H + +namespace Grid { + +template +class CommunicationAvoidingGeneralisedMinimalResidual : public OperatorFunction { + public: + bool ErrorOnNoConverge; // Throw an assert when CAGMRES fails to converge, + // defaults to true + + RealD Tolerance; + + Integer MaxIterations; + Integer RestartLength; + Integer IterationCount; // Number of iterations the CAGMRES took to finish, + // filled in upon completion + + GridStopWatch MatrixTimer; + GridStopWatch LinalgTimer; + GridStopWatch QrTimer; + GridStopWatch CompSolutionTimer; + + Eigen::MatrixXcd H; + + std::vector> y; + std::vector> gamma; + std::vector> c; + std::vector> s; + + CommunicationAvoidingGeneralisedMinimalResidual(RealD tol, + Integer maxit, + Integer restart_length, + bool err_on_no_conv = true) + : Tolerance(tol) + , MaxIterations(maxit) + , RestartLength(restart_length) + , ErrorOnNoConverge(err_on_no_conv) + , H(Eigen::MatrixXcd::Zero(RestartLength, RestartLength + 1)) // sizes taken from DD-αAMG code base + , y(RestartLength + 1, 0.) + , gamma(RestartLength + 1, 0.) + , c(RestartLength + 1, 0.) + , s(RestartLength + 1, 0.) {}; + + void operator()(LinearOperatorBase &LinOp, const Field &src, Field &psi) { + + psi.checkerboard = src.checkerboard; + conformable(psi, src); + + RealD guess = norm2(psi); + assert(std::isnan(guess) == 0); + + RealD cp; + RealD ssq = norm2(src); + RealD rsq = Tolerance * Tolerance * ssq; + + Field r(src._grid); + + std::cout << std::setprecision(4) << std::scientific << std::endl; + std::cout << GridLogIterative << "CommunicationAvoidingGeneralisedMinimalResidual: guess " << guess << std::endl; + std::cout << GridLogIterative << "CommunicationAvoidingGeneralisedMinimalResidual: src " << ssq << std::endl; + + + MatrixTimer.Reset(); + LinalgTimer.Reset(); + QrTimer.Reset(); + CompSolutionTimer.Reset(); + + GridStopWatch SolverTimer; + SolverTimer.Start(); + + IterationCount = 0; + for (int k=0; k &LinOp, const Field &src, Field &psi, RealD rsq) { + + RealD cp = 0; + + Field w(src._grid); + Field r(src._grid); + + std::vector v(RestartLength + 1, src._grid); + + MatrixTimer.Start(); + LinOp.Op(psi, w); + MatrixTimer.Stop(); + + LinalgTimer.Start(); + r = src - w; + + gamma[0] = sqrt(norm2(r)); + + v[0] = (1. / gamma[0]) * r; + LinalgTimer.Stop(); + + for (int i=0; i &LinOp, std::vector &v, Field &w, int iter) { + + MatrixTimer.Start(); + LinOp.Op(v[iter], w); + MatrixTimer.Stop(); + + LinalgTimer.Start(); + for (int i = 0; i <= iter; ++i) { + H(iter, i) = innerProduct(v[i], w); + w = w - H(iter, i) * v[i]; + } + + H(iter, iter + 1) = sqrt(norm2(w)); + v[iter + 1] = (1. / H(iter, iter + 1)) * w; + LinalgTimer.Stop(); + } + + void qrUpdate(int iter) { + + QrTimer.Start(); + for (int i = 0; i < iter ; ++i) { + auto tmp = -s[i] * H(iter, i) + c[i] * H(iter, i + 1); + H(iter, i) = std::conj(c[i]) * H(iter, i) + std::conj(s[i]) * H(iter, i + 1); + H(iter, i + 1) = tmp; + } + + // Compute new Givens Rotation + ComplexD nu = sqrt(std::norm(H(iter, iter)) + std::norm(H(iter, iter + 1))); + c[iter] = H(iter, iter) / nu; + s[iter] = H(iter, iter + 1) / nu; + + // Apply new Givens rotation + H(iter, iter) = nu; + H(iter, iter + 1) = 0.; + + gamma[iter + 1] = -s[iter] * gamma[iter]; + gamma[iter] = std::conj(c[iter]) * gamma[iter]; + QrTimer.Stop(); + } + + void computeSolution(std::vector const &v, Field &psi, int iter) { + + CompSolutionTimer.Start(); + for (int i = iter; i >= 0; i--) { + y[i] = gamma[i]; + for (int k = i + 1; k <= iter; k++) + y[i] = y[i] - H(k, i) * y[k]; + y[i] = y[i] / H(i, i); + } + + if (true) { + for (int i = 0; i <= iter; i++) + psi = psi + v[i] * y[i]; + } + else { + psi = y[0] * v[0]; + for (int i = 1; i <= iter; i++) + psi = psi + v[i] * y[i]; + } + CompSolutionTimer.Stop(); + } +}; +} +#endif diff --git a/tests/solver/Test_wilson_cagmres_unprec.cc b/tests/solver/Test_wilson_cagmres_unprec.cc new file mode 100644 index 00000000..067fc0c1 --- /dev/null +++ b/tests/solver/Test_wilson_cagmres_unprec.cc @@ -0,0 +1,65 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/solver/Test_wilson_cagmres_unprec.cc + +Copyright (C) 2015 + +Author: Daniel Richtmann + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=zero; + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + double volume=1; + for(int mu=0;mu HermOp(Dw); + CommunicationAvoidingGeneralisedMinimalResidual CAGMRES(1.0e-8, 50, 25); + CAGMRES(HermOp,src,result); + + Grid_finalize(); +}