Adding parity flip and gamma 5 option for Block KS starting vector

Grid/algorithms/Algorithms.h

@@ -88,7 +88,7 @@ NAMESPACE_CHECK(multigrid);
 #include <Grid/algorithms/iterative/KrylovSchur.h>
 #include <Grid/algorithms/iterative/BlockKrylovSchur.h>
 #include <Grid/algorithms/iterative/SplitGridBlockKrylovSchur.h>
-//#include <Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h>
+#include <Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h>
 #include <Grid/algorithms/iterative/Gamma5BlockLanczos.h>
 //#include <Grid/algorithms/iterative/Gamma5ScalarLanczos.h>
 #include <Grid/algorithms/iterative/Arnoldi.h>

Grid/algorithms/iterative/BlockKrylovSchur.h

@@ -124,7 +124,17 @@ protected:
 
 public:
   std::vector<Field> evecs;
-  bool doEvalCheck = false;
+  bool doEvalCheck  = false;
+  // When true (and Nblock even), only Nblock/2 starting vectors are required.
+  // The remaining Nblock/2 slots are filled by pairing each supplied vector
+  // with its parity-flipped partner (sign negated on odd checkerboard sites).
+  bool useParityFlip = false;
+  // Same pairing behaviour but the partner is produced by gamma5Func(v).
+  // gamma5Func must be set by the caller before operator() is called.
+  // It is a std::function so that BlockKrylovSchur.h does not need to know
+  // about the Grid QCD Gamma class (which is defined at a later include layer).
+  bool useGamma5     = false;
+  std::function<void(const Field&, Field&)> gamma5Func;
 
   //--------------------------------------------------------------------
   // Constructor
@@ -161,10 +171,17 @@ public:
     Nstop   = _Nstop;
     Nblock  = _Nblock;
 
-    assert((int)v0.size() >= Nblock);
+    {
+      int divisor = 1;
+      if (useParityFlip) divisor *= 2;
+      if (useGamma5)     divisor *= 2;
+      assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor);
+       std::cout << GridLogMessage << "BlockKrylovSchur: divisor = " <<divisor << std::endl;
+    }
     assert(Nm % Nblock == 0);
     assert(Nk % Nblock == 0);
     assert(Nk < Nm);
+    if (useGamma5) assert(gamma5Func && "useGamma5: gamma5Func must be set");
     preRun();  // hook: derived classes add parameter assertions here
 
     int N = Nm;   // total Krylov dimension
@@ -179,8 +196,35 @@ public:
     H = CMat::Zero(N, N);
     B = CMat::Zero(N, Nblock);
 
+    // Build the starting block by expanding each seed with the requested
+    // symmetry partners.  Order: {v, flip(v)} then doubled by {·, g5(·)}
+    // giving {v, flip(v), g5(v), g5(flip(v))} when both flags are set.
+    int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1);
+    std::cout << GridLogMessage << "BlockKrylovSchur: divisor = " <<divisor << std::endl;
     int start = 0;
-    std::vector<Field> startBlock(v0.begin(), v0.begin() + Nblock);
+    std::vector<Field> startBlock;
+    startBlock.reserve(Nblock);
+    for (int i = 0; i < Nblock / divisor; i++) {
+      std::vector<Field> group;
+      group.push_back(v0[i]);
+      if (useParityFlip) {
+        int n = (int)group.size();
+        for (int j = 0; j < n; j++) {
+          Field fp(Grid_);
+          parityFlippedField(group[j], fp);
+          group.push_back(std::move(fp));
+        }
+      }
+      if (useGamma5) {
+        int n = (int)group.size();
+        for (int j = 0; j < n; j++) {
+          Field g5v(Grid_);
+          gamma5Func(group[j], g5v);
+          group.push_back(std::move(g5v));
+        }
+      }
+      for (auto& f : group) startBlock.push_back(std::move(f));
+    }
 
     for (int iter = 0; iter < MaxIter; iter++) {
       std::cout << GridLogMessage << "BlockKrylovSchur: restart iteration " << iter << std::endl;
@@ -253,7 +297,7 @@ public:
       if (Nconv >= Nstop || iter == MaxIter - 1) {
         std::cout << GridLogMessage << "BlockKrylovSchur: done after " << iter
                   << " restarts, " << Nconv << " converged." << std::endl;
-        std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl;
+        std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl;
 
         if (doEvalCheck) {
           Field w(Grid_);
@@ -279,6 +323,21 @@ public:
   CVec                getEvals()         { return evals;         }
   std::vector<RealD>  getRitzEstimates() { return ritzEstimates; }
 
+  // Negate the field on all odd-parity (checkerboard) sites.
+  // Uses pickCheckerboard/setCheckerboard so the SIMD layout is handled
+  // correctly and the negation runs in parallel.
+  static void parityFlippedField(const Field& v, Field& out) {
+    GridCartesian* cgrid = dynamic_cast<GridCartesian*>(v.Grid());
+    assert(cgrid != nullptr);
+    GridRedBlackCartesian rbgrid(cgrid);
+    Field veven(&rbgrid), vodd(&rbgrid);
+    pickCheckerboard(Even, veven, v);
+    pickCheckerboard(Odd,  vodd,  v);
+    vodd = -vodd;
+    setCheckerboard(out, veven);
+    setCheckerboard(out, vodd);
+  }
+
   //--------------------------------------------------------------------
   // Verification: print H and B, check A V = V H + F B^dag explicitly
   //--------------------------------------------------------------------

Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h

@@ -126,7 +126,10 @@ class HarmonicBlockKrylovSchur {
 
 public:
   std::vector<Field>  evecs;
-  bool doEvalCheck = false;
+  bool doEvalCheck  = false;
+  bool useParityFlip = false;
+  bool useGamma5     = false;
+  std::function<void(const Field&, Field&)> gamma5Func;
 
   //--------------------------------------------------------------------
   // Constructor
@@ -153,10 +156,16 @@ public:
     Nstop   = _Nstop;
     Nblock  = _Nblock;
 
-    assert((int)v0.size() >= Nblock);
+    {
+      int divisor = 1;
+      if (useParityFlip) divisor *= 2;
+      if (useGamma5)     divisor *= 2;
+      assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor);
+    }
     assert(Nm % Nblock == 0);
     assert(Nk % Nblock == 0);
     assert(Nk < Nm);
+    if (useGamma5) assert(gamma5Func && "useGamma5: gamma5Func must be set");
 
     int N = Nm;
 
@@ -170,8 +179,31 @@ public:
     H = CMat::Zero(N, N);
     B = CMat::Zero(N, Nblock);
 
+    int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1);
     int start = 0;
-    std::vector<Field> startBlock(v0.begin(), v0.begin() + Nblock);
+    std::vector<Field> startBlock;
+    startBlock.reserve(Nblock);
+    for (int i = 0; i < Nblock / divisor; i++) {
+      std::vector<Field> group;
+      group.push_back(v0[i]);
+      if (useParityFlip) {
+        int n = (int)group.size();
+        for (int j = 0; j < n; j++) {
+          Field fp(Grid_);
+          parityFlippedField(group[j], fp);
+          group.push_back(std::move(fp));
+        }
+      }
+      if (useGamma5) {
+        int n = (int)group.size();
+        for (int j = 0; j < n; j++) {
+          Field g5v(Grid_);
+          gamma5Func(group[j], g5v);
+          group.push_back(std::move(g5v));
+        }
+      }
+      for (auto& f : group) startBlock.push_back(std::move(f));
+    }
 
     for (int iter = 0; iter < MaxIter; iter++) {
       std::cout << GridLogMessage
@@ -280,6 +312,18 @@ public:
   CVec                getEvals()         { return evals;         }
   std::vector<RealD>  getRitzEstimates() { return ritzEstimates; }
 
+  static void parityFlippedField(const Field& v, Field& out) {
+    GridCartesian* cgrid = dynamic_cast<GridCartesian*>(v.Grid());
+    assert(cgrid != nullptr);
+    GridRedBlackCartesian rbgrid(cgrid);
+    Field veven(&rbgrid), vodd(&rbgrid);
+    pickCheckerboard(Even, veven, v);
+    pickCheckerboard(Odd,  vodd,  v);
+    vodd = -vodd;
+    setCheckerboard(out, veven);
+    setCheckerboard(out, vodd);
+  }
+
   //--------------------------------------------------------------------
   // Verification: check A V = V H + F B^dag explicitly
   //--------------------------------------------------------------------

examples/Example_krylov_schur.cc

@@ -334,12 +334,13 @@ int main (int argc, char ** argv)
 
   // Run KrylovSchur and Arnoldi on a Hermitian matrix
     RealD shift=LanParams.shift;
-#if 1
+#if 0
     std::cout << GridLogMessage << "Running Krylov Schur" << std::endl;
     KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
 //    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
     KrySchur.doEvalCheck=true;
-    KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
+//    KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
+    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
     std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
 #else
     std::cout << GridLogMessage << "Running BlockKrylovSchur" << std::endl;
@@ -347,9 +348,16 @@ int main (int argc, char ** argv)
     Nblock=LanParams.Nblock;
     bool if_verify=false;
     if(LanParams.verify) if_verify=true;
-    BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+//    BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
-//    HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall);
+    HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall);
     KrySchur.doEvalCheck=true;
+    KrySchur.useParityFlip=true; std::cout << GridLogMessage << "useParityFlip= " <<KrySchur.useParityFlip<< std::endl;
+    KrySchur.useGamma5=true; std::cout << GridLogMessage << "useGamma5= " <<KrySchur.useGamma5<< std::endl;
+    KrySchur.gamma5Func = [](const FermionField& v, FermionField& out) {
+      Gamma g5(Gamma::Algebra::Gamma5);
+      out = g5 * v;
+    };
+
     KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify);
     std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
 #endif

examples/LanParams.xml

@@ -5,14 +5,14 @@
     <mstep>-0.025</mstep>
     <M5>1.8</M5>
     <Ls>48</Ls>
-    <Nstop>800</Nstop>
+    <Nstop>100</Nstop>
-    <Nk>800</Nk>
+    <Nk>100</Nk>
     <Np>100</Np>
     <ReadEvec>0</ReadEvec>
     <maxIter>1000</maxIter>
     <reorthog>1</reorthog>
     <Nblock>4</Nblock>
-    <verify>1</verify>
+    <verify>0</verify>
     <shift>1.5</shift>
     <resid>1e-8</resid>
     <ChebyLow>1</ChebyLow>