From 841e59e8c00270ea6fdcc9d61b6bd57917f4f129 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 26 Aug 2026 14:22:00 -0400 Subject: [PATCH] GCR coefficient recording and playback --- Grid/algorithms/iterative/GCRCoefficients.h | 59 +++++ Grid/algorithms/multigrid/Smoothers.h | 269 ++++++++++++++++++++ 2 files changed, 328 insertions(+) create mode 100644 Grid/algorithms/iterative/GCRCoefficients.h create mode 100644 Grid/algorithms/multigrid/Smoothers.h diff --git a/Grid/algorithms/iterative/GCRCoefficients.h b/Grid/algorithms/iterative/GCRCoefficients.h new file mode 100644 index 000000000..c09d92747 --- /dev/null +++ b/Grid/algorithms/iterative/GCRCoefficients.h @@ -0,0 +1,59 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/algorithms/iterative/GCRCoefficients.h + + Copyright (C) 2026 + +Author: Peter Boyle + + 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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////////////// +// Recorded GCR coefficients: per-step means over calls of the step length +// a_k and the orthogonalisation coefficients b_kj (already scaled and +// signed as applied: p_{k+1} = r + sum_j b_kj p_{k-j}). +////////////////////////////////////////////////////////////////////////////// +struct GCRCoefficients { + int mmax = 0; + std::vector a_sum; // [k] + std::vector a_n; + std::vector > b_sum; // [k][j] + std::vector > b_n; + void RecordA(int k, ComplexD a){ + if ( (int)a_sum.size() <= k ) { a_sum.resize(k+1,ComplexD(0.0)); a_n.resize(k+1,0); } + a_sum[k] += a; a_n[k]++; + } + void RecordB(int k, const std::vector &b){ + if ( (int)b_sum.size() <= k ) { b_sum.resize(k+1); b_n.resize(k+1); } + if ( b_sum[k].size() < b.size() ) { b_sum[k].resize(b.size(),ComplexD(0.0)); b_n[k].resize(b.size(),0); } + for(int j=0;j<(int)b.size();j++){ b_sum[k][j] += b[j]; b_n[k][j]++; } + } + int Steps(void) const { return a_sum.size(); } + int Calls(void) const { return a_n.size() ? a_n[0] : 0; } + ComplexD A(int k) const { return a_sum[k]/(double)a_n[k]; } + int NB(int k) const { return (k<(int)b_sum.size()) ? b_sum[k].size() : 0; } + ComplexD B(int k,int j) const { return b_sum[k][j]/(double)b_n[k][j]; } + void Report(const std::string &name) const { + std::cout << GridLogMessage << "GCRCoefficients " << name << ": " << Calls() << " calls, " << Steps() << " steps, mmax " << mmax << std::endl; + for(int k=0;k + + 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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once +#include + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////////////// +// Multigrid smoothers as LinearFunction objects. +// +// Gathered here from ~15 copies in tests/solver and tests/debug (HDCR / +// HDCG era), plus the fixed-polynomial smoothers of 2026: +// +// ChebyshevSmoother Chebyshev approx to 1/x on [lo,hi], applied +// through HermOp (the original HDCG smoother on +// a Hermitian shifted operator). +// ChebyshevNonHermitianSmoother same polynomial applied through Op(), for +// a non-Hermitian (near-normal, real-spectrum) +// smoother operator such as shifted PVdagM. +// ChebyshevInverter one Chebyshev-corrected step with residual print. +// MirsSmoother shifted-MdagM CG, HDCG arXiv:1402.2585. +// GCRReplaySmoother replays a GCR's recorded step lengths a_k and +// orthogonalisation coefficients b_kj with NO +// inner products: one matvec per step, zero +// reductions. The "PreconditionerMirsPoly" +// idea of 1402.2585 p.13 applied to GCR. +// +// Recording: PrecGeneralisedConjugateResidualNonHermitian::SetCoefficientRecorder +// accumulates per-step means over calls into a GCRCoefficients; construct a +// GCRReplaySmoother from it. +////////////////////////////////////////////////////////////////////////////// + +inline RealD InverseApproximation(RealD x){ return 1.0/x; } + +////////////////////////////////////////////////////////////////////////////// +// HermOp-based Chebyshev smoother. The second template parameter and the +// 5-argument constructor exist only so the historical call sites +// ChebyshevSmoother S(lo,hi,ord,HermOp,Ddwf); +// compile unchanged; the Matrix argument was never used. +////////////////////////////////////////////////////////////////////////////// +template class ChebyshevSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + typedef LinearOperatorBase FineOperator; + FineOperator & _SmootherOperator; + Chebyshev Cheby; + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) : + _SmootherOperator(SmootherOperator), + Cheby(_lo,_hi,_ord,InverseApproximation) + { + std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator, M &) : + ChebyshevSmoother(_lo,_hi,_ord,SmootherOperator) {}; + void operator() (const Field &in, Field &out) + { + Cheby(_SmootherOperator,in,out); + } +}; + +////////////////////////////////////////////////////////////////////////////// +// Op()-based Chebyshev smoother: x = S(A) r with S the Chebyshev fit to 1/x +// on [lo,hi]. Same three-term recurrence as Chebyshev::operator() +// but through Op, for the non-Hermitian smoother operators of the PVdagM +// multigrid (real coefficients / near-normal, as the recorded GCR +// coefficients show). `order` matvecs, no reductions. +////////////////////////////////////////////////////////////////////////////// +template class ChebyshevNonHermitianSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + LinearOperatorBase &Linop; + RealD lo, hi; int order; + std::vector Coeffs; + ChebyshevNonHermitianSmoother(RealD _lo,RealD _hi,int _order,LinearOperatorBase &Op) + : Linop(Op), lo(_lo), hi(_hi), order(_order) + { + GRID_ASSERT(order>=2); + Coeffs.resize(order); + for(int j=0;j class ChebyshevInverter : public LinearFunction +{ +public: + using LinearFunction::operator(); + typedef LinearOperatorBase FineOperator; + FineOperator & _Operator; + Chebyshev Cheby; + ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) : + _Operator(Operator), + Cheby(_lo,_hi,_ord,InverseApproximation) + { + std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< class MirsSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + typedef LinearOperatorBase FineOperator; + Matrix & SmootherMatrix; + FineOperator & SmootherOperator; + RealD tol; + RealD shift; + int maxit; + MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) : + shift(_shift),tol(_tol),maxit(_maxit), + SmootherOperator(_SmootherOperator), + SmootherMatrix(_SmootherMatrix) + {}; + void operator() (const Field &in, Field &out) + { + ZeroGuesser Guess; + ConjugateGradient CG(tol,maxit,false); + Field src(in.Grid()); + ShiftedMdagMLinearOperator,Field> MdagMOp(SmootherMatrix,shift); + SmootherOperator.AdjOp(in,src); + Guess(src,out); + CG(MdagMOp,src,out); + } +}; + +////////////////////////////////////////////////////////////////////////////// +// Replay of a recorded GCR with a trivial preconditioner: +// p_0 = r_0 ; x_{k+1} = x_k + a_k p_k ; r_{k+1} = r_k - a_k A p_k ; +// p_{k+1} = r_{k+1} + sum_{j class GCRReplaySmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + LinearOperatorBase &Linop; + int mmax, nstep; + std::vector a; + std::vector > b; + GridBase *hist_grid = nullptr; + std::vector p; + GCRReplaySmoother(LinearOperatorBase &Op, const GCRCoefficients &c) : Linop(Op) + { + mmax = c.mmax; GRID_ASSERT(mmax>=1); + nstep = c.Steps(); GRID_ASSERT(nstep>=1); + a.resize(nstep); b.resize(nstep); + for(int k=0;k class SwitchableSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + LinearFunction *current; + std::string label; + SwitchableSmoother(LinearFunction &initial, const std::string &l="initial") : current(&initial), label(l) {} + void Set(LinearFunction &f, const std::string &l) + { + current = &f; label = l; + std::cout << GridLogMessage << " SwitchableSmoother -> " << l << std::endl; + } + void operator() (const Field &in, Field &out) { (*current)(in,out); } +}; + +NAMESPACE_END(Grid);