Blocked->Block. SplitGridBlockKrylovSchur

Grid/algorithms/Algorithms.h

@@ -55,6 +55,7 @@ NAMESPACE_CHECK(approx);
 #include <Grid/algorithms/blas/MomentumProject.h>
 NAMESPACE_CHECK(deflation);
 #include <Grid/algorithms/iterative/ConjugateGradient.h>
+#include <Grid/algorithms/iterative/EigCG.h>
 NAMESPACE_CHECK(ConjGrad);
 #include <Grid/algorithms/iterative/BiCGSTAB.h>
 NAMESPACE_CHECK(BiCGSTAB);
@@ -85,8 +86,9 @@ NAMESPACE_CHECK(multigrid);
 #include <Grid/algorithms/FFT.h>
 
 #include <Grid/algorithms/iterative/KrylovSchur.h>
-#include <Grid/algorithms/iterative/BlockedKrylovSchur.h>
+#include <Grid/algorithms/iterative/BlockKrylovSchur.h>
-#include <Grid/algorithms/iterative/HarmonicBlockedKrylovSchur.h>
+#include <Grid/algorithms/iterative/SplitGridBlockKrylovSchur.h>
+#include <Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h>
 #include <Grid/algorithms/iterative/Arnoldi.h>
 #include <Grid/algorithms/iterative/LanczosBidiagonalization.h>
 #include <Grid/algorithms/iterative/RestartedLanczosBidiagonalization.h>

Grid/algorithms/iterative/BlockKrylovSchur.h

@@ -72,6 +72,7 @@ NAMESPACE_BEGIN(Grid);
 template<class Field>
 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<Field>& 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<Field> W(Nblock, Field(Grid_));
-    for (int t = 0; t < Nblock; t++) {
+    applyBlock(W, kBase);
-      Linop.Op(basis[kBase + t], W[t]);
-    }
 
     // Orthogonalise W against all current basis vectors (full reorthogonalisation)
     // H[i, kBase + t] = <basis[i] | W[t]>

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 <paboyle@ph.ed.ac.uk>
+Author: Chulwoo Jung <chulwoo@bnl.gov>
+
+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<Field> 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 Field>
+class SplitGridBlockKrylovSchur : public BlockKrylovSchur<Field> {
+
+  using Base = BlockKrylovSchur<Field>;
+
+  // Bring protected base members into scope
+  using Base::Nblock;
+  using Base::Grid_;
+  using Base::basis;
+
+  // Split-grid extras
+  LinearOperatorBase<Field>& SLinop;
+  OperatorFunction<Field>&   poly;
+  GridBase*                  SGrid;
+  int                        mrhs;
+
+public:
+
+  SplitGridBlockKrylovSchur(LinearOperatorBase<Field>& _Linop,
+                             LinearOperatorBase<Field>& _SLinop,
+                             OperatorFunction<Field>&   _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<Field>& W, int kBase) override
+  {
+    std::vector<Field> 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

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 <paboyle@ph.ed.ac.uk>
+    Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+    Author: Patrick Oare <patrickoare@gmail.com>
+    Author: Chulwoo Jung <chulwoo@bnl.gov>
 
     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<LatticeFermion> src(Nu,FGrid); 
-  for(int i=0;i<Nu;i++) random(RNG5,src[i]);
+  for(int i=0;i<Nu;i++){
+     random(RNG5,src[i]);
+#if 0
+     LatticeFermion src_e(FrbGrid);
+     pickCheckerboard(Even, src_e, src[i]);
+     src_e=Zero();	
+     setCheckerboard(src[i],src_e);
+     std::cout << GridLogMessage << "src["<<i<<"] even sites zeroed " << std::endl;
+#endif
+  }
 
   if(LanParams.ReadEvec) {
     std::string evecs_file="evec_in";
@@ -336,20 +347,21 @@ int main (int argc, char ** argv)
 //  KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall);
 //  Hacked, really EvalImagSmall
     RealD shift=1.5;
-#if 0
+#if 1
     KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
-    KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
+    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
+//    KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
 #else
     int Nblock=4;
     Nblock=LanParams.Nblock;
     bool if_verify=false;
     if(LanParams.verify) if_verify=true;
-//    KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+    KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
-//    KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,true);
+    KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,true);
-//    BlockedKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+//    BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
 //    KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify);
-    HarmonicBlockedKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall);
+//    HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall);
-    KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true);
+//    KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true);
 #endif
   std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl;