From 88ea24687f731d2c831fbb4d8bb978bb0dce5ad5 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 2 Apr 2026 12:25:55 -0400 Subject: [PATCH] Blocked->Block. SplitGridBlockKrylovSchur --- Grid/algorithms/Algorithms.h | 6 +- Grid/algorithms/iterative/BlockKrylovSchur.h | 33 ++++- .../iterative/SplitGridBlockKrylovSchur.h | 124 ++++++++++++++++++ examples/Example_krylov_schur.cc | 32 +++-- 4 files changed, 179 insertions(+), 16 deletions(-) create mode 100644 Grid/algorithms/iterative/SplitGridBlockKrylovSchur.h diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 6c30beaad..1910a76e9 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -55,6 +55,7 @@ NAMESPACE_CHECK(approx); #include NAMESPACE_CHECK(deflation); #include +#include NAMESPACE_CHECK(ConjGrad); #include NAMESPACE_CHECK(BiCGSTAB); @@ -85,8 +86,9 @@ NAMESPACE_CHECK(multigrid); #include #include -#include -#include +#include +#include +#include #include #include #include diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index 8881bce86..e12604a70 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -72,6 +72,7 @@ NAMESPACE_BEGIN(Grid); template class BlockKrylovSchur { +protected: //-------------------------------------------------------------------- // Types //-------------------------------------------------------------------- @@ -160,6 +161,7 @@ public: assert((int)v0.size() >= Nblock); assert(Nk < Nm); + preRun(); // hook: derived classes add parameter assertions here int N = Nm * Nblock; // total Krylov dimension @@ -403,6 +405,31 @@ public: << "======== end verify ========" << std::endl; } + //-------------------------------------------------------------------- + // Extension hooks for derived classes + //-------------------------------------------------------------------- + + /** + * Apply the operator (or poly-filtered operator) to the current Arnoldi + * block and return the results in W. + * + * Base implementation: W[t] = Linop.Op(basis[kBase + t]). + * Derived classes (e.g. SplitGridBlockKrylovSchur) override this to + * apply poly(A) via Grid_split / Grid_unsplit instead. + */ + virtual void applyBlock(std::vector& W, int kBase) + { + for (int t = 0; t < Nblock; t++) + Linop.Op(basis[kBase + t], W[t]); + } + + /** + * Called once at the start of operator(), after parameters are set. + * Base implementation is a no-op; derived classes override to add + * parameter assertions (e.g. Nblock % mrhs == 0). + */ + virtual void preRun() {} + private: //-------------------------------------------------------------------- @@ -498,11 +525,9 @@ private: int prevN = kBase + Nblock; // number of basis vectors so far after this step int N = Nm * Nblock; - // W[t] = A * basis[kBase + t] + // W[t] = op(basis[kBase + t]) — dispatches to applyBlock() virtual std::vector W(Nblock, Field(Grid_)); - for (int t = 0; t < Nblock; t++) { - Linop.Op(basis[kBase + t], W[t]); - } + applyBlock(W, kBase); // Orthogonalise W against all current basis vectors (full reorthogonalisation) // H[i, kBase + t] = diff --git a/Grid/algorithms/iterative/SplitGridBlockKrylovSchur.h b/Grid/algorithms/iterative/SplitGridBlockKrylovSchur.h new file mode 100644 index 000000000..ae9d14fdc --- /dev/null +++ b/Grid/algorithms/iterative/SplitGridBlockKrylovSchur.h @@ -0,0 +1,124 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/algorithms/iterative/SplitGridBlockKrylovSchur.h + +Copyright (C) 2015 + +Author: Peter Boyle +Author: Chulwoo Jung + +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_SPLIT_BLOCK_KRYLOV_SCHUR_H +#define GRID_SPLIT_BLOCK_KRYLOV_SCHUR_H + +NAMESPACE_BEGIN(Grid); + +/** + * Block Krylov-Schur eigensolver with Split-Grid batched operator application. + * + * Derives from BlockKrylovSchur and overrides only the operator + * application: instead of calling Linop.Op() once per vector, mrhs vectors + * are packed onto a smaller split grid with Grid_split, the polynomial- + * filtered operator poly(SLinop) is applied once to the combined field, and + * the results are unpacked with Grid_unsplit. All other algorithmic logic + * (Arnoldi orthogonalisation, Schur restart, convergence test, etc.) is + * inherited unchanged from the base class. + * + * Constructor extras (beyond BlockKrylovSchur) + * ------------------------------------------- + * SLinop : split-grid linear operator used inside poly + * poly : polynomial filter OperatorFunction applied on the split grid + * SGrid : split-grid GridBase (mrhs fields packed here) + * mrhs : RHS batched per split-grid call; Nblock must be divisible by mrhs + * + * Notes + * ----- + * - The Rayleigh quotient H is built from poly(A), so its Ritz values are + * eigenvalues of poly(A). Convergence is assessed against the full-grid + * Linop (inherited from base), giving true A-eigenvalues. + * - Grid_ (inherited) is the full grid; SGrid is the split grid. + */ +template +class SplitGridBlockKrylovSchur : public BlockKrylovSchur { + + using Base = BlockKrylovSchur; + + // Bring protected base members into scope + using Base::Nblock; + using Base::Grid_; + using Base::basis; + + // Split-grid extras + LinearOperatorBase& SLinop; + OperatorFunction& poly; + GridBase* SGrid; + int mrhs; + +public: + + SplitGridBlockKrylovSchur(LinearOperatorBase& _Linop, + LinearOperatorBase& _SLinop, + OperatorFunction& _poly, + GridBase* _FGrid, + GridBase* _SGrid, + int _mrhs, + RealD _Tolerance, + RitzFilter _rf = EvalReSmall) + : Base(_Linop, _FGrid, _Tolerance, _rf), + SLinop(_SLinop), poly(_poly), SGrid(_SGrid), mrhs(_mrhs) + {} + +protected: + + // Validate mrhs divisibility before the Arnoldi loop starts + void preRun() override + { + assert(Nblock % mrhs == 0 && "Nblock must be divisible by mrhs"); + } + + // Apply poly(A) to basis[kBase .. kBase+Nblock-1] via Grid_split / Grid_unsplit, + // batching mrhs vectors per split-grid call. + void applyBlock(std::vector& W, int kBase) override + { + std::vector in(mrhs, Field(Grid_)); + Field s_in(SGrid); + Field s_out(SGrid); + + int k_start = 0; + while (k_start < Nblock) { + for (int u = 0; u < mrhs; u++) + in[u] = basis[kBase + k_start + u]; + + Grid_split(in, s_in); + poly(SLinop, s_in, s_out); + Grid_unsplit(in, s_out); + + for (int u = 0; u < mrhs; u++) + W[k_start + u] = in[u]; + + k_start += mrhs; + } + } +}; + +NAMESPACE_END(Grid); + +#endif // GRID_SPLIT_BLOCK_KRYLOV_SCHUR_H diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 078588975..c945f852d 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -2,11 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid - Source file: ./tests/Test_padded_cell.cc + Source file: ./examples/Example_krylov_schur.cc Copyright (C) 2023 -Author: Peter Boyle + Author: Peter Boyle + Author: Patrick Oare + Author: Chulwoo Jung 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 @@ -311,7 +313,16 @@ int main (int argc, char ** argv) Nm = Nk + Np; int Nu=16; std::vector src(Nu,FGrid); - for(int i=0;i