Reverted restart to inline in KrylovSchur

Making constructUR to use basisRotate()

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,16 @@ NAMESPACE_CHECK(multigrid);
 #include <Grid/algorithms/FFT.h>
 
 #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/TrueHarmonicBlockKrylovSchur.h>
+#include <Grid/algorithms/iterative/Gamma5BlockLanczos.h>
+//#include <Grid/algorithms/iterative/Gamma5ScalarLanczos.h>
 #include <Grid/algorithms/iterative/Arnoldi.h>
 #include <Grid/algorithms/iterative/LanczosBidiagonalization.h>
 #include <Grid/algorithms/iterative/RestartedLanczosBidiagonalization.h>
+//#include <Grid/algorithms/iterative/GCR.h>
+#include <Grid/algorithms/iterative/MultiSplittingPreconditionedCG.h>
 
 #endif

Grid/algorithms/iterative/BlockKrylovSchur.h

@@ -0,0 +1,878 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/algorithms/iterative/BlockKrylovSchur.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_BLOCKED_KRYLOV_SCHUR_H
+#define GRID_BLOCKED_KRYLOV_SCHUR_H
+
+#include <iomanip>
+#include <numeric>
+
+NAMESPACE_BEGIN(Grid);
+
+/**
+ * Block (block-Arnoldi) restarted Krylov-Schur eigensolver for
+ * general non-Hermitian operators.
+ *
+ * Algorithm
+ * ---------
+ * Uses a block Arnoldi factorisation of block size Nblock:
+ *
+ *   A V_k = V_k H_k + F_k B_k^dag
+ *
+ * where
+ *   V_k   = Nm orthonormal basis vectors (stored flat in basis[])
+ *   H_k   = Nm x Nm upper block-Hessenberg Rayleigh quotient
+ *   F_k   = Nblock residual vectors (the next block beyond V_k)
+ *   B_k   = Nm x Nblock coupling matrix (non-zero only in last Nblock rows)
+ *
+ * Each block Arnoldi step applies A to each of the Nblock vectors in the
+ * current block, orthogonalises against all previous basis vectors, and
+ * reduces the residual block to upper-triangular form via Householder QR
+ * (implemented here as modified Gram-Schmidt within the block).
+ *
+ * The restart is a thick restart via the Schur decomposition of H_k:
+ *   H_k = Q^dag S Q
+ * The leading Nk Schur vectors (chosen by RitzFilter) are retained,
+ * the basis and Rayleigh quotient are truncated, and block Arnoldi continues
+ * from the (Nk/Nblock)-th block.
+ *
+ * Parameters
+ * ----------
+ * Nblock  : block size p
+ * Nm      : total Krylov dimension (must be divisible by Nblock)
+ * Nk      : total vectors to keep after each restart (must be divisible by Nblock, Nk < Nm)
+ * Nstop   : declare convergence when this many eigenpairs have converged
+ * MaxIter : maximum number of outer (restart) iterations
+ * Tolerance : relative convergence tolerance (||r|| < Tolerance * |lambda_max|)
+ */
+template<class Field>
+class BlockKrylovSchur {
+
+protected:
+  //--------------------------------------------------------------------
+  // Types
+  //--------------------------------------------------------------------
+  typedef Eigen::MatrixXcd CMat;
+  typedef Eigen::VectorXcd CVec;
+
+  //--------------------------------------------------------------------
+  // Parameters (set by operator())
+  //--------------------------------------------------------------------
+  int   Nblock;    // block size
+  int   Nm;        // total Krylov dimension (multiple of Nblock)
+  int   Nk;        // total vectors retained after restart (multiple of Nblock)
+  int   Nstop;
+  int   MaxIter;
+  RealD Tolerance;
+
+  //--------------------------------------------------------------------
+  // Internal state
+  //--------------------------------------------------------------------
+  LinearOperatorBase<Field>& Linop;
+  GridBase*                  Grid_;
+  RitzFilter                 ritzFilter;
+
+  // Prefix used in log messages; derived classes override in their constructor
+  // (e.g. "HarmonicBlockKrylovSchur") so the shared code reports the right name.
+  std::string className = "BlockKrylovSchur";
+
+  // Flat storage: basis[s*Nblock + t] is the t-th vector of block s
+  // After construction: basis has Nm entries
+  std::vector<Field> basis;
+
+  // Rayleigh quotient  Nm x Nm
+  CMat H;
+
+  // Residual block: Nblock vectors (the (Nm/Nblock+1)-th block, unnormalised before
+  // QR; normalised and orthogonalised as part of block Arnoldi)
+  std::vector<Field> F;
+
+  // Coupling matrix B: Nm x Nblock.
+  // In exact arithmetic only the last Nblock rows are non-zero:
+  //   B(Nm - Nblock + t, s) = H_{Nm/Nblock+1, Nm/Nblock}(t, s)   (the subdiagonal block)
+  // We keep it as a full matrix for generality after restarts.
+  CMat B;
+
+  RealD beta_k;          // Frobenius norm of the last subdiagonal block
+  RealD rtol;            // absolute tolerance = Tolerance * approxLambdaMax
+
+  // Output
+  CVec               evals;
+  CMat               littleEvecs;   // Nm columns
+  std::vector<RealD> ritzEstimates;
+
+public:
+  std::vector<Field> evecs;
+  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
+  //--------------------------------------------------------------------
+  BlockKrylovSchur(LinearOperatorBase<Field>& _Linop, GridBase* _Grid,
+                     RealD _Tolerance, RitzFilter _rf = EvalReSmall)
+    : Linop(_Linop), Grid_(_Grid), Tolerance(_Tolerance), ritzFilter(_rf),
+      Nblock(-1), Nm(-1), Nk(-1), Nstop(-1), MaxIter(-1),
+      beta_k(0.0), rtol(0.0)
+  {}
+
+  virtual ~BlockKrylovSchur() = default;
+
+  //--------------------------------------------------------------------
+  // Main entry point
+  //--------------------------------------------------------------------
+  /**
+   * Run the blocked Krylov-Schur algorithm.
+   *
+   * Parameters
+   * ----------
+   * v0       : block of Nblock starting vectors (size >= Nblock)
+   * _maxIter : maximum outer (restart) iterations
+   * _Nm      : total Krylov dimension (must be divisible by _Nblock)
+   * _Nk      : total vectors to keep after restart (must be divisible by _Nblock, _Nk < _Nm)
+   * _Nstop   : stop after _Nstop eigenvalues converged
+   * _Nblock  : block size
+   */
+  virtual void operator()(const std::vector<Field>& v0, int _maxIter, int _Nm, int _Nk,
+                          int _Nstop, int _Nblock = 1, bool doubleOrthog = true,
+                          bool doVerify = false)
+  {
+    MaxIter = _maxIter;
+    Nm      = _Nm;
+    Nk      = _Nk;
+    Nstop   = _Nstop;
+    Nblock  = _Nblock;
+
+    {
+      int divisor = 1;
+      if (useParityFlip) divisor *= 2;
+      if (useGamma5)     divisor *= 2;
+      assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor);
+      std::cout << GridLogMessage << className << ": 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
+
+    // Approximate largest eigenvalue for tolerance normalisation
+    RealD approxLambdaMax = approxMaxEval(v0[0]);
+    rtol = Tolerance * approxLambdaMax;
+    std::cout << GridLogMessage << className << ": approx max eval = "
+              << approxLambdaMax << ", rtol = " << rtol << std::endl;
+
+    // Initialise
+    H = CMat::Zero(N, N);
+    B = CMat::Zero(N, Nblock);
+
+    int start = 0;
+    std::vector<Field> startBlock = expandStartBlock(v0);
+
+    for (int iter = 0; iter < MaxIter; iter++) {
+      std::cout << GridLogMessage << className << ": restart iteration " << iter << std::endl;
+
+      // ---- Block Arnoldi: extend from block start to block Nm/Nblock ----
+      blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog);
+
+      // After first full cycle start from block Nk/Nblock
+      start = Nk/Nblock;
+
+      if (doVerify) {
+        std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi";
+        verify(lbl);
+      }
+
+      // ---- Restart rotation: compute Q and the rotated+reordered Rayleigh
+      //      quotient Hnew.  Base class uses the plain Schur decomposition of
+      //      H; derived classes (e.g. HarmonicBlockKrylovSchur) override this
+      //      to target a shift.
+      CMat Q, Hnew;
+      restartRotation(Q, Hnew);
+      CMat Qt = Q.adjoint();
+
+      // Rotate Krylov basis: basis_new[i] = sum_j basis[j] * Qt(j,i)
+      std::vector<Field> basis2;
+      constructUR(basis2, basis, Qt, N);
+      basis = basis2;
+
+      // Update B and H
+      B = Q * B;
+      H = Hnew;
+
+      // ---- Truncate to Nk ----
+      int Nkeep = Nk;
+
+      CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
+      H = CMat::Zero(N, N);
+      H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)) = Htmp;
+
+      std::vector<Field> basisTmp(basis.begin(), basis.begin() + Nkeep);
+      basis = basisTmp;
+
+      CMat Btmp = B(Eigen::seqN(0, Nkeep), Eigen::all);
+      B = CMat::Zero(N, Nblock);
+      B(Eigen::seqN(0, Nkeep), Eigen::all) = Btmp;
+
+      // beta_k = Frobenius norm of the effective coupling
+      beta_k = Btmp.norm();
+      std::cout << GridLogMessage << className << ": beta_k = " << beta_k << std::endl;
+
+      // Restart: the new starting block is F (the residual block from Arnoldi)
+      startBlock = F;
+
+      if (doVerify) {
+        std::string lbl = "iter " + std::to_string(iter) + " after restart+truncation";
+        verify(lbl);
+      }
+
+      // ---- Compute eigensystem of truncated H for convergence check ----
+      CMat Hk = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
+      computeEigensystem(Hk, Nkeep);
+
+      int Nconv = converged(Nkeep);
+      std::cout << GridLogMessage << className << ": converged " << Nconv
+                << " / " << Nstop << std::endl;
+
+      if (Nconv >= Nstop || iter == MaxIter - 1) {
+        std::cout << GridLogMessage << className << ": done after " << iter
+                  << " restarts, " << Nconv << " converged." << std::endl;
+        std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl;
+
+        if (doEvalCheck) {
+          Field w(Grid_);
+          for (int k = 0; k < (int)evecs.size(); k++) {
+            Linop.Op(evecs[k], w);
+            ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+            w -= eval_est * evecs[k];
+            RealD res = std::sqrt(norm2(w));
+            std::cout << GridLogMessage << className << ": evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+
+        return;
+      }
+    }
+  }
+
+  // Accessors
+  std::vector<Field>  getEvecs()         { return evecs;         }
+  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
+  //--------------------------------------------------------------------
+  /**
+   * Checks the block Arnoldi / Krylov-Schur decomposition
+   *
+   *   A V = V H + F B^dag                                      (KS)
+   *
+   * by explicit operator applications.  For each basis vector j:
+   *
+   *   w_j = A basis[j]
+   *
+   * The nBasis × nBasis matrix M of inner products is computed:
+   *
+   *   M[i, j] = <basis[i] | A basis[j]>
+   *
+   * and compared column-by-column against H.  Separately, the nBasis × Nblock
+   * residual coupling matrix R is computed:
+   *
+   *   R[j, t] = <basis[j] | F[t]> * ||F[t]||   (scaled by F-block norms)
+   *
+   * but since F is already normalised, R[j,t] = <basis[j] | F[t]>.
+   *
+   * The KS relation for column j reads:
+   *   w_j = sum_i basis[i] H[i,j]  +  sum_t F[t] B[j,t]*
+   * so the deviation in column j is
+   *   dev_j = w_j - sum_i basis[i] M[i,j]   (should be zero for exact arithmetic)
+   * augmented by the F B^dag term in the last block.
+   *
+   * Prints:
+   *   - H (current Rayleigh quotient, nBasis × nBasis)
+   *   - B (coupling matrix, nBasis × Nblock)
+   *   - M (explicit inner product matrix <V | A V>)
+   *   - max |H[i,j] - M[i,j]|  (should be O(machine epsilon))
+   *   - for each basis column j: || A v_j - V H[:,j] - F B[j,:]^* ||
+   *
+   * Parameters
+   * ----------
+   * label : string printed at the start (e.g. "after restart 2")
+   */
+  void verify(const std::string& label = "")
+  {
+    int nBasis = (int)basis.size();
+    int nF     = (int)F.size();
+
+    if (nBasis == 0) {
+      std::cout << GridLogMessage
+                << className << "::verify [" << label
+                << "]: basis is empty." << std::endl;
+      return;
+    }
+
+    // Full allocated dimension of H (= Nm); may be larger than nBasis after restart+truncation
+    int Nfull = (int)H.rows();
+
+    std::cout << GridLogMessage
+              << "======== " << className << "::verify [" << label << "] ========" << std::endl;
+    std::cout << GridLogMessage
+              << "  nBasis = " << nBasis << "  Nfull = " << Nfull
+              << "  Nblock = " << Nblock << "  nF = " << nF << std::endl;
+
+    // ---- Print H (full Nm x Nm matrix) ----
+    // After restart+truncation nBasis < Nfull; entries outside [0:nBasis,0:nBasis]
+    // should be ~0 — printing the full matrix makes this visible.
+    std::cout << GridLogMessage << "H (" << Nfull << " x " << Nfull << "):" << std::endl;
+    for (int i = 0; i < Nfull; i++) {
+      for (int j = 0; j < Nfull; j++)
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << H(i, j);
+      std::cout << std::endl;
+    }
+
+    // ---- Print B (full Nm x Nblock matrix) ----
+    std::cout << GridLogMessage << "B (" << Nfull << " x " << nF << "):" << std::endl;
+    for (int i = 0; i < Nfull; i++) {
+      for (int t = 0; t < nF; t++)
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << B(i, t);
+      std::cout << std::endl;
+    }
+
+    // ---- Compute M[i,j] = <basis[i] | A basis[j]> ----
+    CMat M = CMat::Zero(nBasis, nBasis);
+    Field w(Grid_);
+    for (int j = 0; j < nBasis; j++) {
+      Linop.Op(basis[j], w);
+      for (int i = 0; i < nBasis; i++)
+        M(i, j) = toStdCmplx(innerProduct(basis[i], w));
+    }
+
+    std::cout << GridLogMessage << "M = <V|AV> (" << nBasis << " x " << nBasis << "):" << std::endl;
+    for (int i = 0; i < nBasis; i++) {
+      for (int j = 0; j < nBasis; j++)
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << M(i, j);
+      std::cout << std::endl;
+    }
+
+    // ---- max |H - M| (for the nBasis x nBasis Rayleigh-quotient block) ----
+    RealD maxHM = 0.0;
+    for (int i = 0; i < nBasis; i++)
+      for (int j = 0; j < nBasis; j++)
+        maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j)));
+    std::cout << GridLogMessage
+              << "  max |H[i,j] - M[i,j]| (i,j < nBasis) = " << maxHM << std::endl;
+
+    // ---- Check entries of H OUTSIDE the nBasis x nBasis block ----
+    // These live in rows i >= nBasis or columns j >= nBasis.
+    // They must be ~0 after restart+truncation (before the coupling rows are
+    // re-filled by the next blockArnoldiIteration call).
+    // During a full Arnoldi run (nBasis == Nfull) this loop is a no-op.
+    RealD maxOutside = 0.0;
+    for (int i = 0; i < Nfull; i++)
+      for (int j = 0; j < Nfull; j++)
+        if (i >= nBasis || j >= nBasis)
+          maxOutside = std::max(maxOutside, std::abs(H(i,j)));
+    if (Nfull > nBasis)
+      std::cout << GridLogMessage
+                << "  max |H[i,j]| outside [0:" << nBasis << ",0:" << nBasis
+                << "] block = " << maxOutside
+                << " (should be ~0 after restart+truncation)" << std::endl;
+
+    // ---- Check orthonormality of basis ----
+    CMat G = CMat::Zero(nBasis, nBasis);
+    for (int i = 0; i < nBasis; i++)
+      for (int j = 0; j < nBasis; j++)
+        G(i, j) = toStdCmplx(innerProduct(basis[i], basis[j]));
+    CMat Gerr = G - CMat::Identity(nBasis, nBasis);
+    std::cout << GridLogMessage
+              << "  max |<V_i|V_j> - delta_ij| = " << Gerr.cwiseAbs().maxCoeff() << std::endl;
+
+    // ---- Per-column residual: || A v_j - V H[:,j] - F B[j,:]^* || ----
+    // For each basis vector j, compute A v_j then subtract V H[:,j] and F B[j,:]^*
+    RealD maxColDev = 0.0;
+    for (int j = 0; j < nBasis; j++) {
+      Linop.Op(basis[j], w);
+
+      // subtract V H[:,j]
+      for (int i = 0; i < nBasis; i++)
+        w -= basis[i] * H(i, j);
+
+      // subtract F B[j,:]^*  (F[t] * conj(B[j,t]))
+      for (int t = 0; t < nF; t++)
+        w -= F[t] * std::conj(B(j, t));
+
+      RealD dev = std::sqrt(norm2(w));
+      std::cout << GridLogMessage
+                << "  || A v[" << j << "] - V H[:,j] - F B[j,:]* || = " << dev << std::endl;
+      maxColDev = std::max(maxColDev, dev);
+    }
+    std::cout << GridLogMessage
+              << "  max column deviation = " << maxColDev << std::endl;
+
+    // ---- Check F block orthogonality against basis ----
+    if (nF > 0) {
+      RealD maxFV = 0.0;
+      for (int t = 0; t < nF; t++)
+        for (int i = 0; i < nBasis; i++) {
+          RealD ip = std::abs(toStdCmplx(innerProduct(basis[i], F[t])));
+          maxFV = std::max(maxFV, ip);
+        }
+      std::cout << GridLogMessage
+                << "  max |<V_i | F_t>| (should be ~0) = " << maxFV << std::endl;
+    }
+
+    std::cout << GridLogMessage
+              << "======== 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() {}
+
+  /**
+   * Compute the restart rotation for the current Rayleigh quotient H.
+   *
+   * Returns Q (unitary) and Hnew = Q H Q^dag reordered so that the wanted
+   * Nk eigenvalues occupy the leading block, ready for truncation.
+   *
+   * Base implementation: the plain Schur decomposition of H, reordered by
+   * ritzFilter.  HarmonicBlockKrylovSchur overrides this to Schur-decompose
+   * the shifted quotient (H - sigma I) so that eigenvalues nearest sigma are
+   * selected.
+   */
+  virtual void restartRotation(CMat& Q, CMat& Hnew)
+  {
+    ComplexSchurDecomposition schur(H, false, ritzFilter);
+    std::cout << GridLogMessage << className << ": Schur decomposed." << std::endl;
+    schur.schurReorder(Nk);
+    std::cout << GridLogMessage << className << ": Schur reordered." << std::endl;
+    Q    = schur.getMatrixQ();
+    Hnew = schur.getMatrixS();
+  }
+
+protected:
+
+  //--------------------------------------------------------------------
+  // Starting-block expansion (parity flip / gamma5 partners)
+  //--------------------------------------------------------------------
+  /**
+   * Expand Nblock/divisor seed vectors into the full Nblock starting block by
+   * pairing each seed with its parity-flipped and/or gamma5 partners.  Order:
+   * {v, flip(v)} then doubled by {., g5(.)} giving {v, flip(v), g5(v),
+   * g5(flip(v))} when both flags are set.  Requires Nblock, Grid_ and the
+   * flags to be set.
+   */
+  std::vector<Field> expandStartBlock(const std::vector<Field>& v0)
+  {
+    int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1);
+    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));
+    }
+    return startBlock;
+  }
+
+protected:
+
+  //--------------------------------------------------------------------
+  // Block Arnoldi iteration
+  //--------------------------------------------------------------------
+  /**
+   * Extends the block Arnoldi factorisation from block index 'start' to
+   * block index 'Nm'.
+   *
+   * On entry (start > 0): basis[0..start*Nblock-1] already set,
+   *   H[0..start*Nblock-1, 0..start*Nblock-1] already set,
+   *   B[start*Nblock-1, :] set (coupling from prior residual block).
+   *   startBlock = the normalised residual block F from the previous cycle.
+   *
+   * On entry (start == 0): initialises everything from startBlock.
+   */
+  void blockArnoldiIteration(std::vector<Field>& startBlock, int endBlock,
+                             int startIdx, bool doubleOrthog)
+  {
+    int N = Nm;
+
+    if (startIdx == 0) {
+      basis.clear();
+      F.clear();
+      H = CMat::Zero(N, N);
+      B = CMat::Zero(N, Nblock);
+
+      // Orthonormalise starting block via modified Gram-Schmidt
+      std::vector<Field> V0 = startBlock;
+      blockOrthonormalise(V0);
+      for (auto& v : V0) basis.push_back(v);
+    } else {
+      // Append residual block (startBlock = F_old) to basis.
+      // The truncated KS relation after restart is:
+      //
+      //   A V_k = V_k S_k + F_old B_old^dag          (*)
+      //
+      // where V_k = basis[0:Nkeep], S_k is stored in H[0:Nkeep,0:Nkeep],
+      // B_old = B[0:Nkeep,:], F_old = startBlock.
+      //
+      // Once F_old is appended as basis[Nkeep:Nkeep+Nblock], (*) becomes
+      // a statement about the extended H matrix:
+      //
+      //   H[Nkeep+t, j] = (B_old^dag)[t,j] = conj(B_old[j,t])
+      //                   for t=0..Nblock-1, j=0..Nkeep-1
+      //
+      // These entries are the "restart coupling rows" that connect the new
+      // block to all retained Schur vectors and must be set before Arnoldi
+      // continues, otherwise A V_k = V_k H[:,0:Nkeep] would be missing the
+      // F_old B_old^dag term for those columns.
+
+      int Nkeep = startIdx * Nblock;
+      for (auto& v : startBlock) basis.push_back(v);
+
+      // Fill restart coupling rows into H
+      for (int t = 0; t < Nblock; t++)
+        for (int j = 0; j < Nkeep; j++)
+          H(Nkeep + t, j) = std::conj(B(j, t));
+
+      // Zero out B for the retained columns now that the coupling is in H
+      for (int j = 0; j < Nkeep; j++)
+        for (int t = 0; t < Nblock; t++)
+          B(j, t) = 0.0;
+    }
+
+    // Main block Arnoldi loop
+    for (int k = startIdx; k < endBlock; k++) {
+      blockArnoldiStep(k, doubleOrthog);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  /**
+   * One block Arnoldi step: extends by one block (Nblock vectors).
+   *
+   * Computes block column k of H and the next basis block V_{k+1}.
+   *
+   * Layout of basis (flat):
+   *   basis[j*Nblock + t]  = t-th vector of j-th block,  j = 0..k
+   *
+   * After this call:
+   *   H[i, k*Nblock : (k+1)*Nblock] filled for i = 0..(k+1)*Nblock - 1
+   *   basis[k*Nblock .. (k+1)*Nblock - 1] normalised (already set on entry)
+   *   F = residual block (to become V_{k+1} after this step if k < Nm-1)
+   *
+   * If k < Nm-1, also:
+   *   H[(k+1)*Nblock : (k+2)*Nblock, k*Nblock : (k+1)*Nblock] = subdiag block (from QR of residual)
+   *   basis extended by Nblock (the normalised residual vectors)
+   */
+  void blockArnoldiStep(int k, bool doubleOrthog)
+  {
+    int kBase  = k * Nblock;          // first flat index of current block
+    int prevN  = kBase + Nblock;      // number of basis vectors so far after this step
+    int N      = Nm;
+
+    // W[t] = op(basis[kBase + t])  — dispatches to applyBlock() virtual
+    std::vector<Field> W(Nblock, Field(Grid_));
+    applyBlock(W, kBase);
+
+    // Orthogonalise W against all current basis vectors (full reorthogonalisation)
+    // H[i, kBase + t] = <basis[i] | W[t]>
+    for (int pass = 0; pass < (doubleOrthog ? 2 : 1); pass++) {
+      for (int i = 0; i < prevN; i++) {
+        for (int t = 0; t < Nblock; t++) {
+          ComplexD coeff = innerProduct(basis[i], W[t]);
+          if (pass == 0)
+            H(i, kBase + t) = toStdCmplx(coeff);
+          else
+            H(i, kBase + t) += toStdCmplx(coeff);
+          W[t] -= coeff * basis[i];
+        }
+      }
+    }
+
+    // Store residual block F
+    F = W;
+
+    if (k == Nm/Nblock - 1) {
+      // Last block: compute coupling matrix B for KS decomp.
+      //
+      // blockQR modifies F in-place (F → Q orthonormal) and returns R
+      // such that  W_orig[t] = sum_s F[s] * R[s,t]  (W_orig = F_after * R).
+      //
+      // The KS relation for column j = kBase+t requires the coefficient of F[s]
+      // to be (B†)[s,j] = conj(B[j,s]).  Matching with R[s,t]:
+      //   conj(B[kBase+t, s]) = R[s,t]  →  B[kBase+t, s] = conj(R[s,t])
+      //
+      // Equivalently the last Nblock rows of B are R^H (Hermitian conjugate of R).
+      // Note: for Nblock=1, R is scalar real positive, so this reduces to B = R. ✓
+      CMat R = blockQR(F);   // F is modified in-place to become Q; returns R
+      for (int t = 0; t < Nblock; t++)
+        for (int s = 0; s < Nblock; s++)
+          B(kBase + t, s) = std::conj(R(s, t));    // B_block = R^H
+
+      beta_k = R.norm();
+      return;
+    }
+
+    // Not last block: QR-decompose residual to get V_{k+1}
+    CMat R = blockQR(F);   // F orthonormalised in-place, R is upper triangular
+
+    // Subdiagonal block of H: H[(k+1)*Nblock : (k+2)*Nblock, kBase : kBase+Nblock] = R
+    int nextBase = (k + 1) * Nblock;
+    for (int i = 0; i < Nblock; i++)
+      for (int j = 0; j < Nblock; j++)
+        H(nextBase + i, kBase + j) = R(i, j);
+
+    // Append normalised residual block to basis
+    for (int t = 0; t < Nblock; t++)
+      basis.push_back(F[t]);
+  }
+
+  //--------------------------------------------------------------------
+  // Block QR via modified Gram-Schmidt within the block
+  //--------------------------------------------------------------------
+  /**
+   * Given a block of Nblock vectors W (not necessarily orthonormal),
+   * orthonormalises them in-place and returns the upper-triangular R
+   * such that W_in = W_out * R.
+   *
+   * Handles (near-)linear dependence by zeroing vectors below threshold.
+   */
+  CMat blockQR(std::vector<Field>& W)
+  {
+    CMat R = CMat::Zero(Nblock, Nblock);
+    const RealD deflThresh = 1e-14;
+
+    for (int j = 0; j < Nblock; j++) {
+      // Orthogonalise W[j] against W[0..j-1]
+      for (int i = 0; i < j; i++) {
+        ComplexD coeff = innerProduct(W[i], W[j]);
+        R(i, j) = toStdCmplx(coeff);
+        W[j] -= coeff * W[i];
+      }
+      RealD nrm = std::sqrt(norm2(W[j]));
+      R(j, j)   = nrm;
+      if (nrm > deflThresh) {
+        W[j] *= (1.0 / nrm);
+      } else {
+        // deflation: zero this vector
+        W[j] = Zero();
+        std::cout << GridLogMessage
+                  << className << ": deflation at block column " << j
+                  << " (norm = " << nrm << ")" << std::endl;
+      }
+    }
+    return R;
+  }
+
+  //--------------------------------------------------------------------
+  // Orthonormalise a block against itself (no prior basis)
+  //--------------------------------------------------------------------
+  void blockOrthonormalise(std::vector<Field>& V)
+  {
+    for (int j = 0; j < (int)V.size(); j++) {
+      for (int i = 0; i < j; i++) {
+        ComplexD c = innerProduct(V[i], V[j]);
+        V[j] -= c * V[i];
+      }
+      RealD nrm = std::sqrt(norm2(V[j]));
+      assert(nrm > 1e-14);
+      V[j] *= (1.0 / nrm);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  // Basis rotation: UR[i] = sum_j U[j] * R(j, i)
+  //--------------------------------------------------------------------
+  void constructUR(std::vector<Field>& UR, std::vector<Field>& U,
+                   CMat& R, int N)
+  {
+    UR.clear();
+    Field tmp(Grid_);
+    for (int i = 0; i < N; i++) {
+      tmp = Zero();
+      for (int j = 0; j < N; j++)
+        tmp += U[j] * R(j, i);
+      UR.push_back(tmp);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  // Eigensystem of the truncated Rayleigh quotient
+  //--------------------------------------------------------------------
+  void computeEigensystem(CMat& Hk, int Nkeep)
+  {
+    Eigen::ComplexEigenSolver<CMat> es;
+    es.compute(Hk);
+
+    // Sort to match schurReorder ordering.
+    int n = es.eigenvalues().size();
+    ComplexComparator cComp(ritzFilter);
+    std::vector<int> idx(n);
+    std::iota(idx.begin(), idx.end(), 0);
+    std::sort(idx.begin(), idx.end(), [&](int a, int b){
+      return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b)));
+    });
+
+    evals.resize(n);
+    littleEvecs.resize(n, n);
+    for (int k = 0; k < n; k++) {
+      evals(k)           = es.eigenvalues()(idx[k]);
+      littleEvecs.col(k) = es.eigenvectors().col(idx[k]);
+    }
+
+    evecs.clear();
+    for (int k = 0; k < Nkeep; k++) {
+      CVec vec = littleEvecs.col(k);
+      Field tmp(Grid_);
+      tmp = Zero();
+      for (int j = 0; j < (int)basis.size() && j < Nkeep; j++)
+        tmp += vec[j] * basis[j];
+      evecs.push_back(tmp);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  // Convergence check
+  //--------------------------------------------------------------------
+  /**
+   * An eigenpair (lambda, y) is converged if the Ritz estimate
+   *   r = || B^dag y ||
+   * satisfies r < rtol.  Here B is the (Nkeep x Nblock) coupling matrix
+   * and y is the little eigenvector (Nkeep-vector) of H.
+   */
+  int converged(int Nkeep)
+  {
+    ritzEstimates.clear();
+    int Nconv = 0;
+
+    CMat Bk = B(Eigen::seqN(0, Nkeep), Eigen::all);  // Nkeep x Nblock
+
+    for (int k = 0; k < Nkeep; k++) {
+      CVec yk   = littleEvecs.col(k);                // Nkeep-vector
+      CVec Bty  = Bk.adjoint() * yk;                // Nblock-vector
+      RealD res = Bty.norm();
+      ritzEstimates.push_back(res);
+      std::cout << GridLogMessage << className << ": Ritz estimate[" << k
+                << "] = " << res << "  eval = " << evals[k] << std::endl;
+      if (res < rtol) Nconv++;
+      else break;
+    }
+    return Nconv;
+  }
+
+  //--------------------------------------------------------------------
+  // Approximate maximum eigenvalue via power iteration
+  //--------------------------------------------------------------------
+  RealD approxMaxEval(const Field& v0, int MAX_ITER = 50)
+  {
+    assert(norm2(v0) > 1e-8);
+    RealD lam = 0.0, denom = std::sqrt(norm2(v0));
+    Field vcur(Grid_), vtmp(Grid_);
+    vcur = v0;
+    for (int i = 0; i < MAX_ITER; i++) {
+      Linop.Op(vcur, vtmp);
+      vcur  = vtmp;
+      RealD num = std::sqrt(norm2(vcur));
+      lam   = num / denom;
+      denom = num;
+    }
+    return lam;
+  }
+
+};
+
+NAMESPACE_END(Grid);
+
+#endif  // GRID_BLOCKED_KRYLOV_SCHUR_H

Grid/algorithms/iterative/Gamma5ScalarLanczos.h

@@ -0,0 +1,230 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    γ5-Scalar Lanczos algorithm for γ5-Hermitian operators (Wilson Dirac).
+
+    Single-vector Lanczos using the γ5-inner product (u,v)_{γ5} = u†γ5v.
+    Builds the Krylov space {q, D_W q, D_W² q, ...} and projects D_W onto it,
+    producing a real tridiagonal matrix T_m whose eigenvalues approximate those
+    of D_W.  Complex conjugate pairs of D_W eigenvalues appear as complex
+    eigenvalues of the real non-symmetric T_m.
+
+    Recurrence (G_k = sign(q_k†γ5 q_k) = ±1):
+        r_k      = D_W q_k − α_k q_k − β_sup[k-1] q_{k-1}
+        β_sub[k] = √|r_k†γ5 r_k|           (lower subdiag T_{k+1,k})
+        β_sup[k] = β_sub[k] G_{k+1}/G_k    (upper superdiag T_{k,k+1})
+        q_{k+1}  = r_k / β_sub[k]
+
+    T_m is real but not symmetric in general; eigenvalues via EigenSolver.
+
+*************************************************************************************/
+#ifndef GRID_GAMMA5_SCALAR_LANCZOS_H
+#define GRID_GAMMA5_SCALAR_LANCZOS_H
+
+#include <functional>
+#include <numeric>
+#include <iomanip>
+
+NAMESPACE_BEGIN(Grid);
+
+template<class Field>
+class Gamma5ScalarLanczos {
+public:
+  using Gamma5Func = std::function<void(const Field&, Field&)>;
+
+private:
+  typedef Eigen::MatrixXcd CMat;
+  typedef Eigen::VectorXcd CVec;
+  typedef Eigen::MatrixXd  RMat;
+
+  LinearOperatorBase<Field>& Linop;
+  GridBase*                  Grid_;
+  Gamma5Func                 applyGamma5;
+  RealD                      Tolerance;
+  int                        nSteps;
+
+  std::vector<Field>  basis;
+  std::vector<RealD>  alpha_;    // diagonal of T_m
+  std::vector<RealD>  beta_sub_; // lower subdiagonal: T_{k+1,k}
+  std::vector<RealD>  beta_sup_; // upper superdiagonal: T_{k,k+1}
+  std::vector<int>    gsign_;    // G_k = ±1
+
+  CVec               evals_;
+  std::vector<Field> evecs_;
+  std::vector<RealD> residuals_;
+
+public:
+  bool doEvalCheck = false;
+
+  Gamma5ScalarLanczos(LinearOperatorBase<Field>& op, GridBase* grid,
+                      Gamma5Func g5, RealD tol = 1e-8)
+    : Linop(op), Grid_(grid), applyGamma5(g5), Tolerance(tol), nSteps(0)
+  {}
+
+  CVec               getEvals()     { return evals_;     }
+  std::vector<Field> getEvecs()     { return evecs_;     }
+  std::vector<RealD> getResiduals() { return residuals_; }
+
+  // v0       : starting vector
+  // maxSteps : maximum Lanczos steps
+  // Nstop    : target Ritz pairs to return (0 = all)
+  // reorthog : full γ5-reorthogonalisation at each step
+  // filter   : eigenvalue selection criterion
+  void operator()(const Field& v0, int maxSteps, int Nstop,
+                  bool reorthog = false, RitzFilter filter = EvalImNormSmall)
+  {
+    basis.clear(); alpha_.clear(); beta_sub_.clear(); beta_sup_.clear(); gsign_.clear();
+    nSteps = 0;
+
+    // Normalise starting vector in γ5-norm: q_0 = v0 / √|v0†γ5 v0|
+    Field q(Grid_); q = v0;
+    Field g5q(Grid_);
+    applyGamma5(q, g5q);
+    RealD omega = std::real(toStdCmplx(innerProduct(g5q, q)));
+    RealD beta0 = std::sqrt(std::abs(omega));
+    assert(beta0 > 1e-14 && "starting vector has zero γ5-norm");
+    q *= (1.0 / beta0);
+    gsign_.push_back(omega >= 0 ? +1 : -1);
+    basis.push_back(q);
+
+    std::cout << GridLogMessage
+              << "Gamma5ScalarLanczos: start  (v0†γ5 v0)=" << omega
+              << "  G[0]=" << gsign_[0] << std::endl;
+
+    for (int k = 0; k < maxSteps; k++) {
+      // ── Apply operator ─────────────────────────────────────────────────
+      Field p(Grid_);
+      Linop.Op(basis[k], p);
+
+      // ── α_k = (q_k†γ5 D_W q_k) / G_k  (real) ─────────────────────────
+      applyGamma5(basis[k], g5q);
+      RealD a = std::real(toStdCmplx(innerProduct(g5q, p))) / RealD(gsign_[k]);
+      alpha_.push_back(a);
+
+      // ── r = D_W q_k − α_k q_k − β_sup[k-1] q_{k-1} ───────────────────
+      Field r(Grid_); r = p;
+      r -= basis[k] * a;
+      if (k > 0)
+        r -= basis[k-1] * beta_sup_[k-1];
+
+      // ── Optional γ5-reorthogonalisation ────────────────────────────────
+      if (reorthog) {
+        for (int i = 0; i <= k; i++) {
+          applyGamma5(basis[i], g5q);
+          ComplexD c = toStdCmplx(innerProduct(g5q, r)) / RealD(gsign_[i]);
+          r -= basis[i] * c;
+        }
+      }
+
+      // ── β_sub[k] = √|r†γ5 r| ───────────────────────────────────────────
+      applyGamma5(r, g5q);
+      omega = std::real(toStdCmplx(innerProduct(g5q, r)));
+      RealD beta = std::sqrt(std::abs(omega));
+
+      std::cout << GridLogMessage
+                << "Gamma5ScalarLanczos: k=" << k
+                << "  alpha=" << a
+                << "  beta=" << beta
+                << "  G[k]=" << gsign_[k]
+                << "  (r†γ5r)=" << omega << std::endl;
+
+      beta_sub_.push_back(beta);
+
+      if (beta < Tolerance) {
+        std::cout << GridLogMessage
+                  << "Gamma5ScalarLanczos: invariant subspace at step " << k << std::endl;
+        nSteps = k + 1;
+        break;
+      }
+
+      int G_next = (omega >= 0) ? +1 : -1;
+      gsign_.push_back(G_next);
+      // β_sup[k] = T_{k,k+1} = β_sub[k] * G_{k+1} / G_k
+      beta_sup_.push_back(beta * RealD(G_next) / RealD(gsign_[k]));
+
+      Field qnew(Grid_); qnew = r; qnew *= (1.0 / beta);
+      basis.push_back(qnew);
+      nSteps = k + 1;
+    }
+
+    if (nSteps == 0) return;
+    computeRitzPairs(nSteps, Nstop, filter);
+  }
+
+private:
+  void computeRitzPairs(int m, int Nstop, RitzFilter filter)
+  {
+    // Assemble real tridiagonal T_m
+    RMat Tm = RMat::Zero(m, m);
+    for (int k = 0; k < m; k++) {
+      Tm(k, k) = alpha_[k];
+      if (k + 1 < m) {
+        Tm(k+1, k) = beta_sub_[k];
+        Tm(k,   k+1) = beta_sup_[k];
+      }
+    }
+
+    std::cout << GridLogMessage
+              << "Gamma5ScalarLanczos: T_m (" << m << "×" << m << "):" << std::endl;
+    for (int i = 0; i < m; i++) {
+      for (int j = 0; j < m; j++)
+        std::cout << "  " << std::setprecision(8) << std::setw(16) << Tm(i,j);
+      std::cout << std::endl;
+    }
+
+    Eigen::EigenSolver<RMat> es(Tm);
+    CVec lambdas = es.eigenvalues();
+    CMat Y       = es.eigenvectors();
+
+    // Sort by filter criterion
+    ComplexComparator cComp(filter);
+    std::vector<int> idx(m);
+    std::iota(idx.begin(), idx.end(), 0);
+    std::sort(idx.begin(), idx.end(), [&](int a, int b){
+      return cComp(toStdCmplx(lambdas(a)), toStdCmplx(lambdas(b)));
+    });
+
+    int Nout = (Nstop > 0 && Nstop < m) ? Nstop : m;
+    evals_.resize(Nout);
+    evecs_.clear();   evecs_.reserve(Nout);
+    residuals_.clear(); residuals_.reserve(Nout);
+
+    // Trailing β for residual estimate: β_sub[m-1] * |y_j[m-1]|
+    RealD beta_trail = (m - 1 < (int)beta_sub_.size()) ? beta_sub_[m-1] : 0.0;
+
+    for (int ji = 0; ji < Nout; ji++) {
+      int  j  = idx[ji];
+      evals_(ji) = lambdas(j);
+      CVec yj = Y.col(j);
+
+      // Ritz vector: u_j = Σ_k q_k y_j[k]
+      Field uj(Grid_); uj = Zero();
+      for (int k = 0; k < m && k < (int)basis.size(); k++)
+        uj += basis[k] * toStdCmplx(yj(k));
+      evecs_.push_back(uj);
+
+      RealD res = beta_trail * std::abs(toStdCmplx(yj(m-1)));
+      if (!std::isfinite(res)) res = std::numeric_limits<RealD>::infinity();
+      residuals_.push_back(res);
+
+      std::cout << GridLogMessage
+                << "  Ritz[" << std::setw(3) << ji << "]"
+                << "  lambda=" << evals_(ji)
+                << "  |res|=" << res;
+
+      if (doEvalCheck) {
+        Field w(Grid_);
+        Linop.Op(uj, w);
+        w -= uj * evals_(ji);
+        RealD actual = std::sqrt(norm2(w));
+        std::cout << "  ||D_W u - λu||=" << actual;
+      }
+      std::cout << std::endl;
+    }
+  }
+};
+
+NAMESPACE_END(Grid);
+
+#endif // GRID_GAMMA5_SCALAR_LANCZOS_H

Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h

@@ -0,0 +1,719 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/algorithms/iterative/HarmonicBlockKrylovSchur.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_HARMONIC_BLOCKED_KRYLOV_SCHUR_H
+#define GRID_HARMONIC_BLOCKED_KRYLOV_SCHUR_H
+
+#include <iomanip>
+#include <numeric>
+
+NAMESPACE_BEGIN(Grid);
+
+/**
+ * Block shift-targeted Krylov-Schur eigensolver.
+ *
+ * Algorithm
+ * ---------
+ * Uses a block Arnoldi factorisation:
+ *
+ *   A V = V H + F B^dag                                         (1)
+ *
+ * with V orthonormal (Nm columns), H the Nm×Nm block upper-Hessenberg
+ * Rayleigh quotient, F the Nblock residual vectors and B the Nm×Nblock
+ * coupling matrix.
+ *
+ * Thick restart
+ * -------------
+ * To target eigenvalues near shift σ, the Schur decomposition is computed
+ * for the shifted Rayleigh quotient:
+ *
+ *   (H - σI) = Q^dag S Q                                        (2)
+ *
+ * Sorting the Schur values of (H - σI) by smallest |S(i,i)| = |λ - σ| and
+ * retaining the leading Nk is equivalent to selecting the Ritz values of H
+ * closest to σ.  Since Q diagonalises H - σI (and hence H itself), the
+ * rotated Rayleigh quotient is exactly upper triangular:
+ *
+ *   H_new = Q H Q^dag = S + σI  (upper triangular)              (3)
+ *
+ * Truncation to Nk is therefore exact: the off-diagonal coupling block
+ * H_new[Nk:, :Nk] = 0 by triangularity.  The Krylov-Schur relation after
+ * restart is exact and block Arnoldi resumes cleanly from F.
+ *
+ * Convergence
+ * -----------
+ * Convergence is declared when || B^H y_k || < Tolerance * approxLambdaMax
+ * for each Ritz pair (λ_k, y_k) of the truncated H.
+ *
+ * Parameters
+ * ----------
+ * shift    : target shift σ (default 0.0); Schur values sorted by |λ - σ|
+ * Nblock   : block size p
+ * Nm       : total Krylov dimension (must be divisible by Nblock)
+ * Nk       : total vectors to retain after each restart (must be divisible by Nblock, Nk < Nm)
+ * Nstop    : stop when this many eigenpairs converge
+ * MaxIter  : maximum outer (restart) iterations
+ * Tolerance: relative convergence tolerance
+ *
+ * Usage
+ * -----
+ *   HarmonicBlockKrylovSchur<Field> hbks(LinOp, Grid, tol, shift, EvalNormSmall);
+ *   std::vector<Field> v0(Nblock, Field(Grid));
+ *   // fill v0 with random starting vectors
+ *   hbks(v0, maxIter, Nm, Nk, Nstop, Nblock);
+ *   auto evals = hbks.getEvals();
+ *   auto evecs = hbks.getEvecs();
+ */
+template<class Field>
+class HarmonicBlockKrylovSchur {
+
+protected:
+  typedef Eigen::MatrixXcd CMat;
+  typedef Eigen::VectorXcd CVec;
+
+  //--------------------------------------------------------------------
+  // Parameters
+  //--------------------------------------------------------------------
+  int    Nblock;
+  int    Nm;        // total Krylov dimension (multiple of Nblock)
+  int    Nk;        // total vectors retained after restart (multiple of Nblock)
+  int    Nstop;
+  int    MaxIter;
+  RealD  Tolerance;
+  ComplexD shift;       // target shift σ
+
+  //--------------------------------------------------------------------
+  // Internal state
+  //--------------------------------------------------------------------
+  LinearOperatorBase<Field>& Linop;
+  GridBase*                  Grid_;
+  RitzFilter                 ritzFilter;
+
+  std::vector<Field> basis;   // Nm*Nblock flat basis
+  CMat               H;       // (Nm*Nblock)² block-Hessenberg Rayleigh quotient
+  std::vector<Field> F;       // Nblock residual vectors
+  CMat               B;       // (Nm*Nblock) × Nblock coupling matrix
+  RealD              beta_k;
+  RealD              rtol;
+
+  CVec               evals;
+  CMat               littleEvecs;
+  std::vector<RealD> ritzEstimates;
+
+public:
+  std::vector<Field>  evecs;
+  bool doEvalCheck  = false;
+  bool useParityFlip = false;
+  bool useGamma5     = false;
+  std::function<void(const Field&, Field&)> gamma5Func;
+
+  //--------------------------------------------------------------------
+  // Constructor
+  //--------------------------------------------------------------------
+  HarmonicBlockKrylovSchur(LinearOperatorBase<Field>& _Linop, GridBase* _Grid,
+                             RealD _Tolerance, ComplexD _shift = 0.0,
+                             RitzFilter _rf = EvalNormSmall)
+    : Linop(_Linop), Grid_(_Grid), Tolerance(_Tolerance), shift(_shift),
+      ritzFilter(_rf),
+      Nblock(-1), Nm(-1), Nk(-1), Nstop(-1), MaxIter(-1),
+      beta_k(0.0), rtol(0.0)
+  {}
+
+  virtual ~HarmonicBlockKrylovSchur() = default;
+
+  //--------------------------------------------------------------------
+  // Main entry point
+  //--------------------------------------------------------------------
+  virtual void operator()(const std::vector<Field>& v0, int _maxIter, int _Nm, int _Nk,
+                          int _Nstop, int _Nblock = 1, bool doubleOrthog = true,
+                          bool doVerify = false)
+  {
+    MaxIter = _maxIter;
+    Nm      = _Nm;
+    Nk      = _Nk;
+    Nstop   = _Nstop;
+    Nblock  = _Nblock;
+
+    {
+      int divisor = 1;
+      if (useParityFlip) divisor *= 2;
+      if (useGamma5)     divisor *= 2;
+      assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor);
+    std::cout << GridLogMessage << "divisor= " << divisor << std::endl;
+
+    }
+    assert(Nm % Nblock == 0);
+    assert(Nk % Nblock == 0);
+    assert(Nk < Nm);
+    if (useGamma5) assert(gamma5Func && "useGamma5: gamma5Func must be set");
+
+    int N = Nm;
+
+    RealD approxLambdaMax = approxMaxEval(v0[0]);
+    rtol = Tolerance * approxLambdaMax;
+    std::cout << GridLogMessage
+              << "HarmonicBlockKrylovSchur: approx max eval = " << approxLambdaMax
+              << ", rtol = " << rtol
+              << ", shift = " << shift << std::endl;
+
+    H = CMat::Zero(N, N);
+    B = CMat::Zero(N, Nblock);
+
+    int start = 0;
+    std::vector<Field> startBlock = expandStartBlock(v0);
+
+    for (int iter = 0; iter < MaxIter; iter++) {
+      std::cout << GridLogMessage
+                << "HarmonicBlockKrylovSchur: restart iteration " << iter << std::endl;
+
+      // ---- Block Arnoldi: extend from block 'start' to block Nm/Nblock ----
+      blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog);
+      start = Nk/Nblock;
+
+      if (doVerify) {
+        std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi";
+        verify(lbl);
+      }
+
+      // ---- Schur decompose (H - σI) to select Schur vectors closest to σ ----
+      // Sorting the Schur values of (H - σI) by |S(i,i)| = |λ - σ| gives the
+      // Ritz values of H nearest the target shift without any matrix inversion.
+      // Because Q diagonalises (H - σI), it also diagonalises H:
+      //   H_new = Q H Q^dag = S + σI  (upper triangular)
+      // Truncation is therefore exact.
+      CMat Hshift = H - shift * CMat::Identity(N, N);
+      ComplexSchurDecomposition schur(Hshift, false, ritzFilter);
+      schur.schurReorder(Nk);
+
+      std::cout << GridLogMessage
+                << "HarmonicBlockKrylovSchur: Ritz values nearest shift (first Nk):" << std::endl;
+      CMat S = schur.getMatrixS();
+      for (int i = 0; i < Nk; i++)
+        std::cout << GridLogMessage << "  [" << i << "] " << S(i, i) + shift << std::endl;
+
+      CMat Q  = schur.getMatrixQ();
+      CMat Qt = Q.adjoint();
+
+      // ---- Rotate Krylov basis ----
+      std::vector<Field> basis2;
+      constructUR(basis2, basis, Qt, N);
+      basis = basis2;
+
+      // ---- Update H and B ----
+      // H_new = S + σI is upper triangular; off-diagonal block H_new[Nk:,:Nk] = 0
+      H = S + shift * CMat::Identity(N, N);
+      B = Q * B;
+
+      // ---- Truncate to Nk (exact: H upper triangular) ----
+      int Nkeep = Nk;
+
+      CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
+      H = CMat::Zero(N, N);
+      H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)) = Htmp;
+
+      std::vector<Field> basisTmp(basis.begin(), basis.begin() + Nkeep);
+      basis = basisTmp;
+
+      CMat Btmp = B(Eigen::seqN(0, Nkeep), Eigen::all);
+      B = CMat::Zero(N, Nblock);
+      B(Eigen::seqN(0, Nkeep), Eigen::all) = Btmp;
+
+      beta_k = Btmp.norm();
+      std::cout << GridLogMessage
+                << "HarmonicBlockKrylovSchur: beta_k = " << beta_k << std::endl;
+
+      // Restart from F (exact: no discarded-basis correction needed)
+      startBlock = F;
+
+      if (doVerify) {
+        std::string lbl = "iter " + std::to_string(iter) + " after restart+truncation";
+        verify(lbl);
+      }
+
+      // ---- Eigensystem of truncated H for convergence ----
+      CMat Hk = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
+      computeEigensystem(Hk, Nkeep);
+
+      int Nconv = converged(Nkeep);
+      std::cout << GridLogMessage
+                << "HarmonicBlockKrylovSchur: converged " << Nconv
+                << " / " << Nstop << std::endl;
+
+      if (Nconv >= Nstop || iter == MaxIter - 1) {
+        std::cout << GridLogMessage
+                  << "HarmonicBlockKrylovSchur: done after " << iter
+                  << " restarts, " << Nconv << " converged." << std::endl;
+        std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl;
+
+        if (doEvalCheck) {
+          Field w(Grid_);
+          for (int k = 0; k < (int)evecs.size(); k++) {
+            Linop.Op(evecs[k], w);
+            ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+            w -= eval_est * evecs[k];
+            RealD res = std::sqrt(norm2(w));
+            std::cout << GridLogMessage << "HarmonicBlockKrylovSchur: evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+
+        return;
+      }
+    }
+  }
+
+  // Accessors
+  std::vector<Field>  getEvecs()         { return evecs;         }
+  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
+  //--------------------------------------------------------------------
+  /**
+   * Checks the block Arnoldi / Krylov-Schur decomposition
+   *
+   *   A V = V H + F B^dag                                      (KS)
+   *
+   * by explicit operator applications.  H here is the standard Rayleigh
+   * quotient, so the KS relation is the same as for BlockKrylovSchur.
+   *
+   * Prints:
+   *   - H (current Rayleigh quotient, nBasis × nBasis)
+   *   - B (coupling matrix, nBasis × Nblock)
+   *   - M (explicit inner product matrix <V | A V>)
+   *   - max |H[i,j] - M[i,j]|  (should be O(machine epsilon))
+   *   - max |<V_i|V_j> - delta_ij|  (orthonormality check)
+   *   - for each basis column j: || A v_j - V H[:,j] - F B[j,:]^* ||
+   *   - max |<V_i | F_t>|  (F orthogonal to basis)
+   */
+  void verify(const std::string& label = "")
+  {
+    int nBasis = (int)basis.size();
+    int nF     = (int)F.size();
+
+    if (nBasis == 0) {
+      std::cout << GridLogMessage
+                << "HarmonicBlockKrylovSchur::verify [" << label
+                << "]: basis is empty." << std::endl;
+      return;
+    }
+
+    std::cout << GridLogMessage
+              << "======== HarmonicBlockKrylovSchur::verify [" << label << "] ========" << std::endl;
+    std::cout << GridLogMessage
+              << "  nBasis = " << nBasis << "  Nblock = " << Nblock
+              << "  nF = " << nF << std::endl;
+
+    // ---- Print H ----
+    std::cout << GridLogMessage << "H (" << nBasis << " x " << nBasis << "):" << std::endl;
+    for (int i = 0; i < nBasis; i++) {
+      for (int j = 0; j < nBasis; j++)
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << H(i, j);
+      std::cout << std::endl;
+    }
+
+    // ---- Print B ----
+    std::cout << GridLogMessage << "B (" << nBasis << " x " << nF << "):" << std::endl;
+    for (int i = 0; i < nBasis; i++) {
+      for (int t = 0; t < nF; t++)
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << B(i, t);
+      std::cout << std::endl;
+    }
+
+    // ---- Compute M[i,j] = <basis[i] | A basis[j]> ----
+    CMat M = CMat::Zero(nBasis, nBasis);
+    Field w(Grid_);
+    for (int j = 0; j < nBasis; j++) {
+      Linop.Op(basis[j], w);
+      for (int i = 0; i < nBasis; i++)
+        M(i, j) = toStdCmplx(innerProduct(basis[i], w));
+    }
+
+    std::cout << GridLogMessage << "M = <V|AV> (" << nBasis << " x " << nBasis << "):" << std::endl;
+    for (int i = 0; i < nBasis; i++) {
+      for (int j = 0; j < nBasis; j++)
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << M(i, j);
+      std::cout << std::endl;
+    }
+
+    // ---- max |H - M| ----
+    RealD maxHM = 0.0;
+    for (int i = 0; i < nBasis; i++)
+      for (int j = 0; j < nBasis; j++)
+        maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j)));
+    std::cout << GridLogMessage
+              << "  max |H[i,j] - M[i,j]| = " << maxHM << std::endl;
+
+    // ---- Check orthonormality of basis ----
+    CMat G = CMat::Zero(nBasis, nBasis);
+    for (int i = 0; i < nBasis; i++)
+      for (int j = 0; j < nBasis; j++)
+        G(i, j) = toStdCmplx(innerProduct(basis[i], basis[j]));
+    CMat Gerr = G - CMat::Identity(nBasis, nBasis);
+    std::cout << GridLogMessage
+              << "  max |<V_i|V_j> - delta_ij| = " << Gerr.cwiseAbs().maxCoeff() << std::endl;
+
+    // ---- Per-column residual: || A v_j - V H[:,j] - F B[j,:]^* || ----
+    RealD maxColDev = 0.0;
+    for (int j = 0; j < nBasis; j++) {
+      Linop.Op(basis[j], w);
+
+      // subtract V H[:,j]
+      for (int i = 0; i < nBasis; i++)
+        w -= basis[i] * H(i, j);
+
+      // subtract F B[j,:]^*  (F[t] * conj(B[j,t]))
+      for (int t = 0; t < nF; t++)
+        w -= F[t] * std::conj(B(j, t));
+
+      RealD dev = std::sqrt(norm2(w));
+      std::cout << GridLogMessage
+                << "  || A v[" << j << "] - V H[:,j] - F B[j,:]* || = " << dev << std::endl;
+      maxColDev = std::max(maxColDev, dev);
+    }
+    std::cout << GridLogMessage
+              << "  max column deviation = " << maxColDev << std::endl;
+
+    // ---- Check F block orthogonality against basis ----
+    if (nF > 0) {
+      RealD maxFV = 0.0;
+      for (int t = 0; t < nF; t++)
+        for (int i = 0; i < nBasis; i++) {
+          RealD ip = std::abs(toStdCmplx(innerProduct(basis[i], F[t])));
+          maxFV = std::max(maxFV, ip);
+        }
+      std::cout << GridLogMessage
+                << "  max |<V_i | F_t>| (should be ~0) = " << maxFV << std::endl;
+    }
+
+    std::cout << GridLogMessage
+              << "======== end verify ========" << std::endl;
+  }
+
+protected:
+
+  //--------------------------------------------------------------------
+  // Starting-block expansion (parity flip / gamma5 partners)
+  //--------------------------------------------------------------------
+  // Expands Nblock/divisor seed vectors into the full Nblock starting
+  // block, pairing each seed with its parity-flipped and/or gamma5
+  // partners.  Requires Nblock, Grid_ and the flags to be set.
+  std::vector<Field> expandStartBlock(const std::vector<Field>& v0)
+  {
+    int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1);
+    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));
+    }
+    return startBlock;
+  }
+
+  //--------------------------------------------------------------------
+  // Block Arnoldi iteration
+  //--------------------------------------------------------------------
+  void blockArnoldiIteration(std::vector<Field>& startBlock, int endBlock,
+                             int startIdx, bool doubleOrthog)
+  {
+    int N = Nm;
+
+    if (startIdx == 0) {
+      basis.clear();
+      F.clear();
+      H = CMat::Zero(N, N);
+      B = CMat::Zero(N, Nblock);
+
+      std::vector<Field> V0 = startBlock;
+      blockOrthonormalise(V0);
+      for (auto& v : V0) basis.push_back(v);
+    } else {
+      // Append the new starting block to the retained basis.
+      //
+      // The exact truncated KS relation is  A V_k = V_k H_k + F B_k^dag,
+      // so the coupling rows are  H[Nkeep+t, j] = conj(B_k[j,t]).
+      // Since H_new = S + σI is upper triangular, the off-diagonal block
+      // H_new[Nkeep:, :Nkeep] = 0 and the restart from F is exact.
+      int Nkeep = startIdx * Nblock;
+      for (auto& v : startBlock) basis.push_back(v);
+
+      // Fill restart coupling rows into H (exact: H_new is upper triangular,
+      // so B encodes the only non-zero coupling to the new block).
+      for (int t = 0; t < Nblock; t++)
+        for (int j = 0; j < Nkeep; j++)
+          H(Nkeep + t, j) = std::conj(B(j, t));
+
+      // Zero out B for the retained columns now that the coupling is in H
+      for (int j = 0; j < Nkeep; j++)
+        for (int t = 0; t < Nblock; t++)
+          B(j, t) = 0.0;
+    }
+
+    for (int k = startIdx; k < endBlock; k++)
+      blockArnoldiStep(k, doubleOrthog);
+  }
+
+  //--------------------------------------------------------------------
+  // One block Arnoldi step
+  //--------------------------------------------------------------------
+  void blockArnoldiStep(int k, bool doubleOrthog)
+  {
+    int kBase  = k * Nblock;
+    int prevN  = kBase + Nblock;
+    int N      = Nm;
+
+    std::vector<Field> W(Nblock, Field(Grid_));
+    for (int t = 0; t < Nblock; t++)
+      Linop.Op(basis[kBase + t], W[t]);
+
+    // Full reorthogonalisation against all current basis vectors
+    for (int pass = 0; pass < (doubleOrthog ? 2 : 1); pass++) {
+      for (int i = 0; i < prevN; i++) {
+        for (int t = 0; t < Nblock; t++) {
+          ComplexD coeff = innerProduct(basis[i], W[t]);
+          if (pass == 0)
+            H(i, kBase + t) = toStdCmplx(coeff);
+          else
+            H(i, kBase + t) += toStdCmplx(coeff);
+          W[t] -= coeff * basis[i];
+        }
+      }
+    }
+
+    F = W;
+
+    if (k == Nm/Nblock - 1) {
+      // Last block: record coupling in B as R^H (Hermitian conjugate of QR factor)
+      // KS relation requires B[kBase+t, s] = conj(R[s,t])
+      CMat R = blockQR(F);
+      for (int t = 0; t < Nblock; t++)
+        for (int s = 0; s < Nblock; s++)
+          B(kBase + t, s) = std::conj(R(s, t));    // B_block = R^H
+      beta_k = R.norm();
+      return;
+    }
+
+    // Not last: QR the residual, extend basis
+    CMat R = blockQR(F);
+
+    int nextBase = (k + 1) * Nblock;
+    for (int i = 0; i < Nblock; i++)
+      for (int j = 0; j < Nblock; j++)
+        H(nextBase + i, kBase + j) = R(i, j);
+
+    for (int t = 0; t < Nblock; t++)
+      basis.push_back(F[t]);
+
+  }
+
+  //--------------------------------------------------------------------
+  // Block QR (modified Gram-Schmidt within the block)
+  //--------------------------------------------------------------------
+  CMat blockQR(std::vector<Field>& W)
+  {
+    CMat R = CMat::Zero(Nblock, Nblock);
+    const RealD deflThresh = 1e-14;
+
+    for (int j = 0; j < Nblock; j++) {
+      for (int i = 0; i < j; i++) {
+        ComplexD coeff = innerProduct(W[i], W[j]);
+        R(i, j) = toStdCmplx(coeff);
+        W[j] -= coeff * W[i];
+      }
+      RealD nrm = std::sqrt(norm2(W[j]));
+      R(j, j)   = nrm;
+      if (nrm > deflThresh) {
+        W[j] *= (1.0 / nrm);
+      } else {
+        W[j] = Zero();
+        std::cout << GridLogMessage
+                  << "HarmonicBlockKrylovSchur: deflation at block column " << j
+                  << " (norm = " << nrm << ")" << std::endl;
+      }
+    }
+    return R;
+  }
+
+  //--------------------------------------------------------------------
+  // Orthonormalise starting block
+  //--------------------------------------------------------------------
+  void blockOrthonormalise(std::vector<Field>& V)
+  {
+    for (int j = 0; j < (int)V.size(); j++) {
+      for (int i = 0; i < j; i++) {
+        ComplexD c = innerProduct(V[i], V[j]);
+        V[j] -= c * V[i];
+      }
+      RealD nrm = std::sqrt(norm2(V[j]));
+      assert(nrm > 1e-14);
+      V[j] *= (1.0 / nrm);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  // Basis rotation: UR[i] = sum_j U[j] * R(j,i)
+  //--------------------------------------------------------------------
+  void constructUR(std::vector<Field>& UR, std::vector<Field>& U,
+                   CMat& R, int N)
+  {
+    UR.clear();
+    Field tmp(Grid_);
+    for (int i = 0; i < N; i++) {
+      tmp = Zero();
+      for (int j = 0; j < N; j++)
+        tmp += U[j] * R(j, i);
+      UR.push_back(tmp);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  // Eigensystem of the truncated H (not Hhat)
+  //--------------------------------------------------------------------
+  /**
+   * Eigenvalues of H_k are the standard Ritz values in the retained
+   * subspace.  After convergence has been declared via harmonic estimates,
+   * the final reported eigenvalues and vectors come from H_k (not Hhat_k),
+   * since H_k contains the true projected operator.
+   */
+  void computeEigensystem(CMat& Hk, int Nkeep)
+  {
+    Eigen::ComplexEigenSolver<CMat> es;
+    es.compute(Hk);
+
+    // Sort to match schurReorder ordering.
+    int n = es.eigenvalues().size();
+    ComplexComparator cComp(ritzFilter);
+    std::vector<int> idx(n);
+    std::iota(idx.begin(), idx.end(), 0);
+    std::sort(idx.begin(), idx.end(), [&](int a, int b){
+      return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b)));
+    });
+
+    evals.resize(n);
+    littleEvecs.resize(n, n);
+    for (int k = 0; k < n; k++) {
+      evals(k)           = es.eigenvalues()(idx[k]);
+      littleEvecs.col(k) = es.eigenvectors().col(idx[k]);
+    }
+
+    evecs.clear();
+    for (int k = 0; k < Nkeep; k++) {
+      CVec vec = littleEvecs.col(k);
+      Field tmp(Grid_);
+      tmp = Zero();
+      for (int j = 0; j < (int)basis.size() && j < Nkeep; j++)
+        tmp += vec[j] * basis[j];
+      evecs.push_back(tmp);
+    }
+  }
+
+  //--------------------------------------------------------------------
+  // Convergence check
+  //--------------------------------------------------------------------
+  /**
+   * Ritz estimate for eigenpair k: || B^H y_k ||
+   * where y_k is the k-th eigenvector of the truncated H.
+   * The same bound applies whether using Ritz or harmonic Ritz restart.
+   */
+  int converged(int Nkeep)
+  {
+    ritzEstimates.clear();
+    int Nconv = 0;
+
+    CMat Bk = B(Eigen::seqN(0, Nkeep), Eigen::all);
+
+    for (int k = 0; k < Nkeep; k++) {
+      CVec yk   = littleEvecs.col(k);
+      CVec Bty  = Bk.adjoint() * yk;
+      RealD res = Bty.norm();
+      ritzEstimates.push_back(res);
+      std::cout << GridLogMessage
+                << "HarmonicBlockKrylovSchur: Ritz estimate[" << k
+                << "] = " << res << "  eval = " << evals[k] << std::endl;
+      if (res < rtol) Nconv++;
+      else break;
+    }
+    return Nconv;
+  }
+
+  //--------------------------------------------------------------------
+  // Approximate maximum eigenvalue (power iteration)
+  //--------------------------------------------------------------------
+  RealD approxMaxEval(const Field& v0, int MAX_ITER = 50)
+  {
+    assert(norm2(v0) > 1e-8);
+    RealD lam = 0.0, denom = std::sqrt(norm2(v0));
+    Field vcur(Grid_), vtmp(Grid_);
+    vcur = v0;
+    for (int i = 0; i < MAX_ITER; i++) {
+      Linop.Op(vcur, vtmp);
+      vcur  = vtmp;
+      RealD num = std::sqrt(norm2(vcur));
+      lam   = num / denom;
+      denom = num;
+    }
+    return lam;
+  }
+
+};
+
+NAMESPACE_END(Grid);
+
+#endif  // GRID_HARMONIC_BLOCKED_KRYLOV_SCHUR_H

Grid/algorithms/iterative/KrylovSchur.h

@@ -9,6 +9,7 @@ Copyright (C) 2015
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: paboyle <paboyle@ph.ed.ac.uk>
 Author: Patrick Oare <poare@bnl.gov>
+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
@@ -42,6 +43,19 @@ using std::conj;
 inline const ComplexD& toStdCmplx(const ComplexD& c) { return c; }
 #endif
 
+/**
+ * Small dense coefficient matrix of ComplexD (thrust::complex on GPU builds)
+ * for passing to basisRotate, whose device kernel cannot multiply Grid tensors
+ * by std::complex. Must be namespace-scope: nvcc rejects function-local types
+ * as template arguments to functions containing extended device lambdas.
+ */
+struct KSCoeffMat {
+  std::vector<ComplexD> a; int n;
+  KSCoeffMat(int _n) : a(_n*_n), n(_n) {}
+  ComplexD&       operator()(int i, int j)       { return a[i*n+j]; }
+  const ComplexD& operator()(int i, int j) const { return a[i*n+j]; }
+};
+
 /**
  * Options for which Ritz values to keep in implicit restart. TODO move this and utilities into a new file
  */
@@ -105,8 +119,8 @@ struct ComplexComparator
   bool operator()(std::complex<double> z1, std::complex<double> z2) { 
 	  RealD tmp1=std::abs(std::imag(z1));
 	  RealD tmp2=std::abs(std::imag(z2));
-	  if ( std::abs(std::real(z1)) >4.) tmp1 += 100.;
-	  if ( std::abs(std::real(z2)) >4.) tmp2 += 100.;
+	  if ( std::abs(std::real(z1)) >2.) tmp1 += 100.;
+	  if ( std::abs(std::real(z2)) >2.) tmp2 += 100.;
     switch (RF) {
       case EvalNormSmall:
         return std::abs(z1) < std::abs(z2);
@@ -336,12 +350,12 @@ class KrylovSchur {
 
     RitzFilter ritzFilter;                  // how to sort evals
 
-  public:       
-    RealD *shift;
+  public:
+    bool doEvalCheck = false;
 
     KrylovSchur(LinearOperatorBase<Field> &_Linop, GridBase *_Grid, RealD _Tolerance, RitzFilter filter = EvalReSmall)
       : Linop(_Linop), Grid(_Grid), Tolerance(_Tolerance), ritzFilter(filter), u(_Grid), MaxIter(-1), Nm(-1), Nk(-1), Nstop (-1),
-        evals (0), ritzEstimates (), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0),shift(NULL)
+        evals (0), ritzEstimates (), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0)
     {
       u = Zero();
     };
@@ -358,34 +372,26 @@ class KrylovSchur {
     std::vector<Field>  getEvecs()             { return evecs;         }
 
     /**
-     * Runs the Krylov-Schur loop.
+     * Runs the non-harmonic Krylov-Schur loop.
      *   - Runs an Arnoldi step to generate the Rayleigh quotient and Krylov basis.
-     *   - Schur decompose the Rayleigh quotient. 
-     *   - Permutes the Rayleigh quotient according to the eigenvalues. 
-     *   - Truncate the Krylov-Schur expansion. 
+     *   - Schur decompose the Rayleigh quotient.
+     *   - Permutes the Rayleigh quotient according to the eigenvalues.
+     *   - Truncate the Krylov-Schur expansion.
      */
-    void operator()(const Field& v0, int _maxIter, int _Nm, int _Nk, int _Nstop, RealD *_shift=NULL, bool doubleOrthog = true) {
-
-//      RealD shift_=1.;
-//      shift = &shift_;
-      if(_shift) shift = _shift;
+    void operator()(const Field& v0, int _maxIter, int _Nm, int _Nk, int _Nstop, bool doubleOrthog = true) {
 
       MaxIter = _maxIter;
       Nm = _Nm; Nk = _Nk;
       Nstop = _Nstop;
-      
+
       ssq = norm2(v0);
       RealD approxLambdaMax = approxMaxEval(v0);
       rtol = Tolerance * approxLambdaMax;
       std::cout << GridLogMessage << "Approximate max eigenvalue: " << approxLambdaMax << std::endl;
-      // rtol = Tolerance;
 
       b = Eigen::VectorXcd::Zero(Nm);       // start as e_{k+1}
       b(Nm-1) = 1.0;
 
-      // basis = new std::vector<Field> (Nm, Grid);
-      // evecs.reserve();
-
       int start = 0;
       Field startVec = v0;
       littleEvecs = Eigen::MatrixXcd::Zero(Nm, Nm);
@@ -394,134 +400,13 @@ class KrylovSchur {
 
         // Perform Arnoldi steps to compute Krylov basis and Rayleigh quotient (Hess)
         arnoldiIteration(startVec, Nm, start, doubleOrthog);
-        startVec = u;        // original code
+        startVec = u;
         start = Nk;
-	 
+
         std::cout << GridLogDebug << "b after Arnoldi " << b << std::endl;
 
-        // checkKSDecomposition();
-
-        RealD gamma;
-        Field uhat(Grid);
-        Eigen::MatrixXcd Btilde;
-	std::vector<Field> basis2_s;
-        Eigen::VectorXcd b_s;
-#if 1
-if (shift){
-
-if(0){
-      Field w(Grid);
-
-      ComplexD coeff,coeff2;
-      for (int j = 0; j < Nm; j++) {
-        Linop.Op(basis[j], w);
-        for (int k = 0; k < Nm; k++) {
-          coeff2 = innerProduct(basis[k], basis[j]);
-          coeff = innerProduct(basis[k], w);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
-        std::cout << GridLogMessage << " Rayleigh "<<k<<" "<<j<<" "<<Rayleigh (k,j)<<" "<<coeff << " <k|j> = " << coeff2 << std::endl;
-        }
-        coeff = innerProduct(basis[j], u);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
-        std::cout << GridLogMessage << " u "<<j<<" "<<coeff << std::endl;
-
-      }
-}
-	Eigen::MatrixXcd temp = Rayleigh;
-        for (int m=0;m<Nm;m++) temp(m,m) -= *shift;
-        Eigen::MatrixXcd RayleighS = temp.inverse(); // (B-tI)^-1
-        Eigen::MatrixXcd temp2;
-
-//        temp2	= RayleighS*temp;
-//        std::cout << GridLogDebug  << "Shift inverse check: shift= "<<*shift<<" "<< temp2 <<std::endl;
-
-        temp2=RayleighS.adjoint(); //(B-tI)^-1*
-        Eigen::VectorXcd g  = temp2*b; //g = (B-tI)^-1* * b
-        Btilde= Rayleigh + g*(b.adjoint());
-
-        Field utilde(Grid);
-        utilde  = u;
-        for (int j = 0; j<Nm; j++){
-          utilde -= basis[j]*g(j);
-        }
-
-	ComplexSchurDecomposition schurS (Btilde, false, ritzFilter);
-	std::cout << GridLogMessage << "Shifted Schur eigenvalues shift = "<<*shift  << std::endl;
-        schurS.schurReorder(Nk);
-
-	Eigen::MatrixXcd Q_s = schurS.getMatrixQ();
-        Eigen::MatrixXcd Qt_s = Q_s.adjoint();                           // TODO should Q be real?
-#if 0
-	std::cout << GridLogMessage << "Q_s" << Q_s <<std::endl;
-	std::cout << GridLogMessage << "Qt_s" << Qt_s <<std::endl;
-	Eigen::MatrixXcd temp4= Q_s*Qt_s;
-	std::cout << GridLogMessage << "Q_s*Qt_s" << temp4 << std::endl;
-	temp4 = Btilde;
-	std::cout << GridLogMessage << "Btilde" << temp4<< std::endl;
-	temp4 = temp4*Qt_s;
-	std::cout << GridLogMessage << "Btilde*Qt_s" << temp4<< std::endl;
-	temp4 = Q_s*temp4;
-	std::cout << GridLogMessage << "Q_s*Btilde*Qt_s" << temp4<< std::endl;
-#endif
-
-        Eigen::MatrixXcd S_s = schurS.getMatrixS();
-
-        Btilde=schurS.getMatrixS();
-        b_s= b;
-        b_s=Q_s*b; // Q is Qt in SlepC, b_s=bhat
-
-        constructUR(basis2_s, basis, Qt_s, Nm);
-
-
-	Eigen::MatrixXcd RayTmp_s = Btilde(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk));
-        Btilde = RayTmp_s;
-
-        std::vector<Field> basisTmp_s = std::vector<Field> (basis2_s.begin(), basis2_s.begin() + Nk);
-        basis2_s = basisTmp_s;
-
-        Eigen::VectorXcd btmp_s = b_s.head(Nk);
-        b_s = btmp_s;
-	
-        Eigen::VectorXcd ghat = g;
-	ghat = -Q_s * g;
-
-        Eigen::VectorXcd gtmp_s = ghat.head(Nk);
-        ghat = gtmp_s;
-
-        uhat  = utilde;
-        for (int j = 0; j<Nk; j++){
-          uhat -= basis2_s[j]*ghat(j);
-        }
-
-        gamma = std::sqrt(norm2(uhat));
-	uhat *= 1.0/gamma;
-        std::cout << GridLogMessage << " gamma "<<gamma << std::endl;
-
-	Btilde  += ghat*(b_s.adjoint());
-	b_s *=gamma;
-
-// Eq.(44)
-    if(0){
-          Field w(Grid);
-    
-          ComplexD coeff,coeff2;
-          for (int j = 0; j < Nk; j++) {
-            Linop.Op(basis2_s[j], w);
-            for (int k = 0; k < Nk; k++) {
-              coeff2 = innerProduct(basis2_s[k], basis2_s[j]);
-              coeff = innerProduct(basis2_s[k], w);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
-            std::cout << GridLogMessage << " Btilde "<<k<<" "<<j<<" "<<Btilde(k,j)<<" "<<coeff << " <k|j> = " << coeff2 << std::endl;
-            }
-            coeff = innerProduct(basis2_s[j], uhat);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
-            coeff2 = innerProduct(uhat,w);
-            std::cout << GridLogMessage << " uhat "<<j<<" "<<coeff << " w "<< coeff2 << " b " << b_s (j) << " ghat "<<ghat(j)<< std::endl;
-    
-          }
-    }
-} 
-#endif
-
-
-if (!shift){
-        // Perform a Schur decomposition on Rayleigh
+        // --- Restart: Schur-decompose Rayleigh, reorder by ritzFilter, rotate
+        //     the basis into the Schur vectors, and truncate to the leading Nk.
         ComplexSchurDecomposition schur (Rayleigh, false, ritzFilter);
         std::cout << GridLogDebug << "Schur decomp holds? " << schur.checkDecomposition() << std::endl;
 
@@ -532,8 +417,6 @@ if (!shift){
 
         Eigen::MatrixXcd Q = schur.getMatrixQ();
         Qt = Q.adjoint();                           // TODO should Q be real?
-        Eigen::MatrixXcd S = schur.getMatrixS();
-        // std::cout << GridLogMessage << "Schur decomp holds after reorder? " << schur.checkDecomposition() << std::endl;
 
         std::cout << GridLogMessage << "*** ROTATING TO SCHUR BASIS *** " << std::endl;
 
@@ -541,78 +424,167 @@ if (!shift){
         Rayleigh = schur.getMatrixS();
         b = Q * b;            // b^\dag = b^\dag * Q^\dag <==> b = Q*b
 
-        std::vector<Field> basis2; 
-        constructUR(basis2, basis, Qt, Nm);
-        basis = basis2;
-if(0){
-      Field w(Grid);
-
-      ComplexD coeff,coeff2;
-      for (int j = 0; j < Nm; j++) {
-        Linop.Op(basis[j], w);
-        for (int k = 0; k < Nm; k++) {
-          coeff2 = innerProduct(basis[k], basis[j]);
-          coeff = innerProduct(basis[k], w);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
-        std::cout << GridLogMessage << " Stilde "<<k<<" "<<j<<" "<<Rayleigh(k,j)<<" "<<coeff << " <k|j> = " << coeff2 << std::endl;
-        }
-        coeff = innerProduct(basis[j], u);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
-        std::cout << GridLogMessage << " u"<<j<<" "<<coeff << std::endl;
-
-      }
-}
-}
-
+        std::vector<Field> basis2;
+        constructUR(basis2, basis, Qt, Nm,Nk);
 
         std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** " << std::endl;
-if (!shift){
         std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl;
 
         Eigen::MatrixXcd RayTmp = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk));
         Rayleigh = RayTmp;
 
-        std::vector<Field> basisTmp = std::vector<Field> (basis.begin(), basis.begin() + Nk);
-        basis = basisTmp;
+        basis = std::vector<Field> (basis2.begin(), basis2.begin() + Nk);
 
-        Eigen::VectorXcd btmp = b.head(Nk);
-        b = btmp;
+        b = b.head(Nk).eval();
 
         std::cout << GridLogDebug << "Rayleigh after truncation: " << std::endl << Rayleigh << std::endl;
 
-
         checkKSDecomposition();
 
         // Compute eigensystem of Rayleigh. Note the eigenvectors correspond to the sorted eigenvalues.
         computeEigensystem(Rayleigh);
         std::cout << GridLogMessage << "Eigenvalues (first Nk sorted): " << std::endl << evals << std::endl;
-}
 
-	if(shift){
-		Rayleigh = Btilde;
-		basis= basis2_s;
-		b  = b_s;
-		beta_k = gamma;
-		u= uhat;
+        if (checkConvergedAndReport(i)) return;
+      }
+    }
+
+    /**
+     * Runs the harmonic (shifted) Krylov-Schur loop: extracts Ritz values of a
+     * shift-augmented Rayleigh quotient so that eigenvalues near `*_shift` are
+     * reordered to the top instead of the extremal ones.
+     */
+    void operator()(const Field& v0, int _maxIter, int _Nm, int _Nk, int _Nstop, RealD *_shift, bool doubleOrthog = true) {
+
+      assert(_shift && "harmonic KrylovSchur: shift must be non-null");
+      RealD shiftVal = *_shift;
+
+      MaxIter = _maxIter;
+      Nm = _Nm; Nk = _Nk;
+      Nstop = _Nstop;
+
+      ssq = norm2(v0);
+      RealD approxLambdaMax = approxMaxEval(v0);
+      rtol = Tolerance * approxLambdaMax;
+      std::cout << GridLogMessage << "Approximate max eigenvalue: " << approxLambdaMax << std::endl;
+
+      b = Eigen::VectorXcd::Zero(Nm);       // start as e_{k+1}
+      b(Nm-1) = 1.0;
+
+      int start = 0;
+      Field startVec = v0;
+      littleEvecs = Eigen::MatrixXcd::Zero(Nm, Nm);
+      for (int i = 0; i < MaxIter; i++) {
+        std::cout << GridLogMessage << "Restart Iteration " << i << std::endl;
+
+        // Perform Arnoldi steps to compute Krylov basis and Rayleigh quotient (Hess)
+        arnoldiIteration(startVec, Nm, start, doubleOrthog);
+        startVec = u;
+        start = Nk;
+
+        std::cout << GridLogDebug << "b after Arnoldi " << b << std::endl;
+
+        // --- Harmonic restart: Schur-decompose the shift-augmented Rayleigh
+        //     quotient so Ritz values close to shiftVal are reordered to the
+        //     top, then rotate and truncate as in the non-harmonic case.
+        Eigen::MatrixXcd temp = Rayleigh;
+        for (int m=0;m<Nm;m++) temp(m,m) -= shiftVal;
+        Eigen::MatrixXcd RayleighS = temp.inverse(); // (B-tI)^-1
+        Eigen::MatrixXcd temp2 = RayleighS.adjoint(); //(B-tI)^-1*
+        Eigen::VectorXcd g  = temp2*b; //g = (B-tI)^-1* * b
+        Eigen::MatrixXcd Btilde = Rayleigh + g*(b.adjoint());
+
+        Field utilde(Grid);
+        utilde  = u;
+        for (int j = 0; j<Nm; j++){
+          utilde -= basis[j]*g(j);
+        }
+
+        ComplexSchurDecomposition schurS (Btilde, false, ritzFilter);
+        std::cout << GridLogMessage << "Shifted Schur eigenvalues shift = "<<shiftVal  << std::endl;
+        schurS.schurReorder(Nk);
+
+        Eigen::MatrixXcd Q_s = schurS.getMatrixQ();
+        Eigen::MatrixXcd Qt_s = Q_s.adjoint();                           // TODO should Q be real?
+
+        Btilde=schurS.getMatrixS();
+        Eigen::VectorXcd b_s = Q_s*b; // Q is Qt in SlepC, b_s=bhat
+
+        std::vector<Field> basis2_s;
+        constructUR(basis2_s, basis, Qt_s, Nm,Nk);
+
+        Eigen::MatrixXcd RayTmp_s = Btilde(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk));
+        Btilde = RayTmp_s;
+
+        basis2_s = std::vector<Field> (basis2_s.begin(), basis2_s.begin() + Nk);
+
+        b_s = b_s.head(Nk).eval();
+
+        Eigen::VectorXcd ghat = -Q_s * g;
+        ghat = ghat.head(Nk).eval();
+
+        Field uhat(Grid);
+        uhat  = utilde;
+        for (int j = 0; j<Nk; j++){
+          uhat -= basis2_s[j]*ghat(j);
+        }
+
+        RealD gamma = std::sqrt(norm2(uhat));
+        uhat *= 1.0/gamma;
+        std::cout << GridLogMessage << " gamma "<<gamma << std::endl;
+
+        Btilde  += ghat*(b_s.adjoint());
+        b_s *=gamma;
+
+        Rayleigh = Btilde;
+        basis= basis2_s;
+        b  = b_s;
+        beta_k = gamma;
+        u= uhat;
+
         checkKSDecomposition();
         computeEigensystem(Rayleigh);
         std::cout << GridLogMessage << "Eigenvalues (first Nk sorted): " << std::endl << evals << std::endl;
-	}
 
-        // check convergence and return if needed. 
-        int Nconv = converged();
-        std::cout << GridLogMessage << "Number of evecs converged: " << Nconv << std::endl;
-        if (Nconv >= Nstop || i == MaxIter - 1) {
-          std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " 
-                        << i << "." << std::endl;
-          // basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm);      // Think this might have been the issue
-          std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl;
-
-          // writeEigensystem(path);
-
-          return;
-        }
+        if (checkConvergedAndReport(i)) return;
       }
     }
 
+  private:
+
+    /**
+     * Shared post-restart convergence check. On convergence (or the final
+     * iteration) reports eigenvalues and optionally the explicit residual
+     * check. Returns true if operator() should return.
+     */
+    bool checkConvergedAndReport(int i) {
+      int Nconv = converged();
+      std::cout << GridLogMessage << "Number of evecs converged: " << Nconv << std::endl;
+      if (Nconv >= Nstop || i == MaxIter - 1) {
+        std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration "
+                      << i << "." << std::endl;
+        std::cout << GridLogMessage << "Eigenvalues: " << std::endl << evals << std::endl;
+
+        if (doEvalCheck) {
+          Field w(Grid);
+          for (int k = 0; k < (int)evecs.size(); k++) {
+            Linop.Op(evecs[k], w);
+            ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+            w -= eval_est * evecs[k];
+            RealD res = std::sqrt(norm2(w));
+            std::cout << GridLogMessage << "KrylovSchur: evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+        return true;
+      }
+      return false;
+    }
+
+  public:
+
     /**
      * Constructs the Arnoldi basis for the Krylov space K_n(D, src). (TODO make private)
      * 
@@ -723,35 +695,36 @@ if (!shift){
     {
       std::cout << GridLogMessage << "Computing eigenvalues." << std::endl;
 
-      // evals = S.diagonal();
-      int n = evals.size();       // should be regular Nm
-
       evecs.clear();
-      // evecs.assign(n, Field(Grid));
-
-      // TODO: is there a faster way to get the eigenvectors of a triangular matrix?
-      // Rayleigh.triangularView<Eigen::Upper> tri;
 
       Eigen::ComplexEigenSolver<Eigen::MatrixXcd> es;
-      // es.compute(Rayleigh);
       es.compute(S);
-      evals = es.eigenvalues();
-      littleEvecs = es.eigenvectors();
 
-      // std::cout << GridLogDebug << "Little evecs: " << littleEvecs << std::endl;
-      // std::cout << "Rayleigh diag: " << S.diagonal() << std::endl;
-      // std::cout << "Rayleigh evals: " << evals << std::endl;
+      // Sort eigenvalues/evecs to match the schurReorder ordering.
+      int n = es.eigenvalues().size();
+      ComplexComparator cComp(ritzFilter);
+      std::vector<int> idx(n);
+      std::iota(idx.begin(), idx.end(), 0);
+      std::sort(idx.begin(), idx.end(), [&](int a, int b){
+        return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b)));
+      });
+
+      evals.resize(n);
+      littleEvecs.resize(n, n);
+      for (int k = 0; k < n; k++) {
+        evals(k)          = es.eigenvalues()(idx[k]);
+        littleEvecs.col(k) = es.eigenvectors().col(idx[k]);
+      }
 
       // Convert evecs to lattice fields
       for (int k = 0; k < n; k++) {
         Eigen::VectorXcd vec = littleEvecs.col(k);
         Field tmp (basis[0].Grid());
         tmp = Zero();
-        for (int j = 0; j < basis.size(); j++) {
+        for (int j = 0; j < (int)basis.size(); j++) {
           tmp = tmp + vec[j] * basis[j];
         }
         evecs.push_back(tmp);
-        // evecs[k] = tmp;
       }
     }
 
@@ -816,12 +789,12 @@ if (!shift){
         Eigen::VectorXcd evec_k = littleEvecs.col(k);
         RealD ritzEstimate = std::abs(b.dot(evec_k));           // b^\dagger s
         ritzEstimates.push_back(ritzEstimate);
-        // ritzEstimates[k] = ritzEstimate;
         std::cout << GridLogMessage << "Ritz estimate for evec " << k << " = " << ritzEstimate << std::endl;
         if (ritzEstimate < rtol) {
           Nconv++;
+        } else {
+          break;
         }
-
       }
       // Check that Ritz estimate is explicitly || D (Uy) - lambda (Uy) ||
 //       checkRitzEstimate();
@@ -844,7 +817,7 @@ if (!shift){
       // std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl;
 
       std::vector<Field> rotated = basis;
-      constructUR(rotated, basis, Rayleigh, k);             // manually rotate
+      constructUR(rotated, basis, Rayleigh, k,k);             // manually rotate
       // Eigen::MatrixXcd Rt = Rayleigh.adjoint();
       // basisRotate(rotated, Rt, 0, k, 0, k, k);           // UR
 
@@ -934,31 +907,24 @@ if (!shift){
     }
 
     /**
-     * Given a vector of fields U (equivalently, a LxN matrix, where L is the number of degrees of 
-     * freedom on the lattice field) and an NxN matrix R, forms the product UR. 
-     * 
-     * Note that I believe this is equivalent to basisRotate(U, R.adjoint(), 0, N, 0, N, N), but I'm 
-     * not 100% sure (this will be slower and unoptimized though).
+     * Given a vector of fields U (equivalently, a LxN matrix, where L is the number of degrees of
+     * freedom on the lattice field) and an NxN matrix R, forms the product UR.
+     * Only the first N2 output columns are computed; columns [N2, N) are Zero.
+     *
+     * Implemented with the accelerator basisRotate kernel, which computes
+     * basis[j] <- sum_k Qt(j,k) basis[k] in place, so we copy U into UR and
+     * pass R^T. The coefficients are staged as ComplexD (thrust::complex on
+     * GPU builds) since std::complex has no device arithmetic with Grid tensors.
      */
-    void constructUR(std::vector<Field>& UR, std::vector<Field> &U, Eigen::MatrixXcd& R, int N) {
-      Field tmp (Grid);
+    void constructUR(std::vector<Field>& UR, std::vector<Field> &U, Eigen::MatrixXcd& R, int N, int N2) {
+      KSCoeffMat Rt(N);
+      for (int i = 0; i < N; i++)
+        for (int j = 0; j < N; j++)
+          Rt(i, j) = ComplexD(R(j, i).real(), R(j, i).imag());
 
-      UR.clear();
-      // UR.resize(N);
-
-      std::cout << GridLogDebug << "R to rotate by (should be Rayleigh): " << R << std::endl;
-
-      for (int i = 0; i < N; i++) {
-        tmp = Zero();
-        for (int j = 0; j < N; j++) {
-          std::cout << GridLogDebug << "Adding R("<<j<<", "<<i<<") = " << R(j, i) << " to rotated" << std::endl;
-          std::cout << GridLogDebug << "Norm of U[j] is " << norm2(U[j]) << " to rotated" << std::endl;
-          tmp = tmp + U[j] * R(j, i);
-        }
-        std::cout << GridLogDebug << "rotated norm at i = " << i << " is: " << norm2(tmp) << std::endl;
-        UR.push_back(tmp);
-        // UR[i] = tmp;
-      }
+      UR = U;
+      basisRotate(UR, Rt, 0, N2, 0, N, N);
+      for (int i = N2; i < N; i++) UR[i] = Zero();
       return;
     }
 

Grid/algorithms/iterative/RestartedLanczosBidiagonalization.h

@@ -153,7 +153,7 @@ public:
       // Step 1: extend from pStart steps to Nm steps
       // ----------------------------------------------------------------
       extendBasis(pStart, Nm, betaRestart);
-      verify();
+//      verify();
 
       // ----------------------------------------------------------------
       // Step 2: SVD of the Nm x Nm B matrix.
@@ -208,7 +208,7 @@ public:
       // Step 4: implicit restart — compress to Nk steps
       // ----------------------------------------------------------------
       implicitRestart(Nm, Nk, sigma, X, Y, order, betaK, betaRestart);
-      verify();
+//      verify();
 
       // Lucky breakdown: exact invariant subspace found; convergence is exact.
       // B_p^+ = diag(alpha[0..Nk-1]); extract directly from restart basis.
@@ -360,7 +360,7 @@ private:
     if (restart_col >= 0 && restart_col < m && (int)fvec.size() > 0) {
       for (int j = 0; j < restart_col && j < (int)fvec.size(); ++j){
         B(j, restart_col) = fvec[j];
-        std::cout << GridLogMessage << "buildFullB: B  " <<j<<" "<<restart_col<<B(j, restart_col);
+        std::cout << GridLogDebug << "buildFullB: B  " <<j<<" "<<restart_col<<B(j, restart_col)<<std::endl;
       }
     }
     return B;
@@ -464,20 +464,150 @@ private:
     }
   }
 
+public:
   // ------------------------------------------------------------------
-  // Full reorthogonalization of vec against the vectors in basis.
-  // Subtracts projections only — does NOT normalize.
+  // Block reorthogonalization helpers.
+  // Declared public because CUDA extended lambdas cannot live inside
+  // private/protected member functions.
+  //
+  // batchInnerProducts: computes c[j] = <basis[j], vec> for all j
+  //   in a single GPU pass (one accelerator_barrier instead of n).
+  //   Queues n pairs of (per-site kernel, reduceKernel) to computeStream
+  //   without intermediate CPU syncs, then syncs once at the end.
+  //
+  // batchUpdate: computes vec -= sum_j c[j]*basis[j] in one GPU kernel.
+  //
+  // reorthogonalize: two-pass Classical Gram-Schmidt (CGS2) using the
+  //   two helpers above.  Each pass costs 2 GPU syncs (1 IP + 1 update)
+  //   instead of 2n syncs per pass in the old sequential MGS.
+  // ------------------------------------------------------------------
+
+  void batchInnerProducts(const Field &vec,
+                          const std::vector<Field> &basis,
+                          std::vector<ComplexD> &c)
+  {
+    int n = (int)basis.size();
+    c.resize(n);
+    if (n == 0) return;
+
+    typedef typename Field::vector_object         vobj;
+    typedef decltype(innerProduct(vobj(), vobj())) inner_t;
+    typedef decltype(basis[0].View(AcceleratorRead)) View;
+
+    GridBase *grid = vec.Grid();
+    uint64_t oSites = grid->oSites();
+    uint64_t nsimd  = grid->Nsimd();
+
+    // all_ip[j * oSites + ss] = per-site inner product of basis[j] and vec at site ss.
+    // Layout: n contiguous blocks of oSites each.
+    deviceVector<inner_t> all_ip((uint64_t)n * oSites);
+    inner_t *all_ip_p = &all_ip[0];
+
+    hostVector<View>   h_basis_v(n);
+    deviceVector<View> d_basis_v(n);
+    for (int j = 0; j < n; ++j) {
+      h_basis_v[j] = basis[j].View(AcceleratorRead);
+      acceleratorPut(d_basis_v[j], h_basis_v[j]);
+    }
+    View *basis_vp = &d_basis_v[0];
+
+    // Queue n per-site kernels to the accelerator stream — no intermediate barriers.
+    autoView(vec_v, vec, AcceleratorRead);
+    for (int j = 0; j < n; ++j) {
+      int      jj      = j;
+      uint64_t oSites_ = oSites;
+      accelerator_for(ss, oSites, nsimd, {
+        auto x = coalescedRead(basis_vp[jj][ss]);
+        auto y = coalescedRead(vec_v[ss]);
+        coalescedWrite(all_ip_p[jj * oSites_ + ss], innerProduct(x, y));
+      });
+    }
+
+    // ONE sync after all n kernels
+    accelerator_barrier();
+
+    // Copy all per-site results to host
+    hostVector<inner_t> all_ip_h((uint64_t)n * oSites);
+    acceleratorCopyFromDevice(all_ip_p, &all_ip_h[0], (uint64_t)n * oSites * sizeof(inner_t));
+
+    // Reduce on host: sum over oSites, then collapse SIMD lanes via Reduce(TensorRemove(...))
+    // TensorRemove strips the iSinglet tensor wrapper to expose the SIMD scalar type.
+    // Reduce sums all nsimd lanes and returns a plain scalar (RealD or ComplexD).
+    std::vector<ComplexD> raw(n);
+    for (int j = 0; j < n; ++j) {
+      inner_t sum = Zero();
+      for (uint64_t ss = 0; ss < oSites; ++ss)
+        sum += all_ip_h[(uint64_t)j * oSites + ss];
+      raw[j] = ComplexD(Reduce(TensorRemove(sum)));
+    }
+    grid->GlobalSumVector(&raw[0], n);
+    for (int j = 0; j < n; ++j) c[j] = raw[j];
+
+    for (int j = 0; j < n; ++j) h_basis_v[j].ViewClose();
+  }
+
+  void batchUpdate(Field &vec,
+                   const std::vector<Field> &basis,
+                   const std::vector<ComplexD> &c)
+  {
+    int n = (int)basis.size();
+    if (n == 0) return;
+
+    typedef typename Field::vector_object vobj;
+    typedef decltype(basis[0].View(AcceleratorRead)) View;
+
+    GridBase *grid = vec.Grid();
+    uint64_t oSites = grid->oSites();
+    uint64_t nsimd  = grid->Nsimd();
+
+    // Split complex coefficients into real/imag double arrays on device.
+    // Using doubles avoids potential ComplexD-device-code compatibility issues.
+    hostVector<double>   h_re(n), h_im(n);
+    deviceVector<double> d_re(n), d_im(n);
+    for (int k = 0; k < n; ++k) {
+      h_re[k] = c[k].real();
+      h_im[k] = c[k].imag();
+    }
+    acceleratorCopyToDevice(&h_re[0], &d_re[0], n * sizeof(double));
+    acceleratorCopyToDevice(&h_im[0], &d_im[0], n * sizeof(double));
+    double *re_p = &d_re[0];
+    double *im_p = &d_im[0];
+
+    // Basis views
+    hostVector<View>   h_basis_v(n);
+    deviceVector<View> d_basis_v(n);
+    for (int k = 0; k < n; ++k) {
+      h_basis_v[k] = basis[k].View(AcceleratorRead);
+      acceleratorPut(d_basis_v[k], h_basis_v[k]);
+    }
+    View *basis_vp = &d_basis_v[0];
+
+    // Single kernel: vec[ss] -= sum_k (re[k] + i*im[k]) * basis[k][ss]
+    autoView(vec_v, vec, AcceleratorWrite);
+    accelerator_for(ss, oSites, nsimd, {
+      auto v = coalescedRead(vec_v[ss]);
+      for (int k = 0; k < n; ++k) {
+        auto b = coalescedRead(basis_vp[k][ss]);
+        v = v - re_p[k] * b - timesI(im_p[k] * b);
+      }
+      coalescedWrite(vec_v[ss], v);
+    });
+
+    for (int k = 0; k < n; ++k) h_basis_v[k].ViewClose();
+  }
+
+  // ------------------------------------------------------------------
+  // Full reorthogonalization using two-pass Classical Gram-Schmidt (CGS2).
+  // Each pass calls batchInnerProducts (1 GPU sync) + batchUpdate (1 sync),
+  // replacing the old 2n GPU syncs per pass from sequential MGS.
   // ------------------------------------------------------------------
   void reorthogonalize(Field &vec, const std::vector<Field> &basis)
   {
-    for (int j = 0; j < (int)basis.size(); ++j) {
-      ComplexD ip = innerProduct(basis[j], vec);
-      vec = vec - ip * basis[j];
-    }
-    // Second pass for numerical stability
-    for (int j = 0; j < (int)basis.size(); ++j) {
-      ComplexD ip = innerProduct(basis[j], vec);
-      vec = vec - ip * basis[j];
+    if (basis.empty()) return;
+    std::vector<ComplexD> c;
+    for (int pass = 0; pass < 2; ++pass) {
+      batchInnerProducts(vec, basis, c);
+      batchUpdate(vec, basis, c);
     }
   }
 

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

Grid/algorithms/iterative/TrueHarmonicBlockKrylovSchur.h

@@ -0,0 +1,366 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/algorithms/iterative/TrueHarmonicBlockKrylovSchur.h
+
+Copyright (C) 2015
+
+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_TRUE_HARMONIC_BLOCKED_KRYLOV_SCHUR_H
+#define GRID_TRUE_HARMONIC_BLOCKED_KRYLOV_SCHUR_H
+
+#include <Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h>
+
+NAMESPACE_BEGIN(Grid);
+
+/**
+ * True harmonic-Ritz block Krylov-Schur eigensolver.
+ *
+ * Difference from the base class
+ * ------------------------------
+ * HarmonicBlockKrylovSchur is a *shift-targeted* Krylov-Schur: it Schur-
+ * decomposes (H - sigma I) and sorts the standard Ritz values by |lambda -
+ * sigma|.  For interior eigenvalues of Hermitian indefinite operators the
+ * standard Ritz values can be spurious ("ghosts"), so the sort can select
+ * junk.  This variant performs genuine harmonic Ritz extraction instead.
+ *
+ * Harmonic extraction
+ * -------------------
+ * From the exact block Arnoldi relation
+ *
+ *   A V = V H + F B^dag                                            (1)
+ *
+ * harmonic Ritz pairs (theta, g) w.r.t. shift sigma satisfy the Petrov-
+ * Galerkin condition (A - sigma)Vg - (theta - sigma)Vg  perp  (A - sigma)V,
+ * which reduces (with Hs = H - sigma I) to the small eigenproblem
+ *
+ *   Hhat g = (theta - sigma) g,   Hhat = Hs + Hs^{-H} B B^dag      (2)
+ *
+ * Exact thick restart
+ * -------------------
+ * From (2), H g = theta g - Hs^{-H} B (B^dag g), so the residual of every
+ * harmonic pair lies in the Nblock-dimensional column space of the single
+ * block
+ *
+ *   W = F - V Hs^{-H} B :
+ *   A (V g) - theta (V g) = W (B^dag g)                            (3)
+ *
+ * Hence {V g_1 .. g_Nk} + cols(W) admits an EXACT truncated Krylov-Schur
+ * relation  A Y = Y H_k + What Bc^dag  with dense H_k, and block Arnoldi
+ * resumes cleanly from What.  Krylov-Schur (Stewart, SIMAX 23(3), 2001)
+ * does not require H to be triangular after restart -- only that the
+ * decomposition holds exactly, which (3) guarantees.
+ *
+ * This is the block analogue of the GMRES-DR restart (Morgan, SISC 24(1),
+ * 2002; block version: Morgan, Appl. Numer. Math. 54, 2005).  The exact
+ * composition "block harmonic Krylov-Schur" is a derivation from these
+ * ingredients rather than a published algorithm; the relation (1) is
+ * checked explicitly at every restart when doVerify=true, and identity (3)
+ * was validated to machine precision in a standalone Eigen prototype.
+ *
+ * Caveats
+ * -------
+ * - If sigma is (numerically) an eigenvalue of H, Hs^{-H} is singular;
+ *   an assert fires -- perturb the shift slightly.
+ * - Reported eigenvalues come from the eigensystem of the retained H_k
+ *   (the true projected Rayleigh quotient Y^dag A Y), as in the base
+ *   class, so convergence tests and reporting are inherited unchanged.
+ *
+ * Usage: identical to HarmonicBlockKrylovSchur.
+ */
+template<class Field>
+class TrueHarmonicBlockKrylovSchur : public HarmonicBlockKrylovSchur<Field> {
+
+  typedef HarmonicBlockKrylovSchur<Field> Base;
+  typedef typename Base::CMat CMat;
+  typedef typename Base::CVec CVec;
+
+  // Members of the dependent base class need explicit scope.
+  using Base::Nblock;
+  using Base::Nm;
+  using Base::Nk;
+  using Base::Nstop;
+  using Base::MaxIter;
+  using Base::Tolerance;
+  using Base::shift;
+  using Base::Linop;
+  using Base::Grid_;
+  using Base::ritzFilter;
+  using Base::basis;
+  using Base::H;
+  using Base::F;
+  using Base::B;
+  using Base::beta_k;
+  using Base::rtol;
+  using Base::evals;
+  using Base::littleEvecs;
+  using Base::ritzEstimates;
+
+  using Base::expandStartBlock;
+  using Base::blockArnoldiIteration;
+  using Base::blockQR;
+  using Base::computeEigensystem;
+  using Base::converged;
+  using Base::approxMaxEval;
+
+public:
+  using Base::evecs;
+  using Base::doEvalCheck;
+  using Base::useParityFlip;
+  using Base::useGamma5;
+  using Base::gamma5Func;
+  using Base::verify;
+
+  TrueHarmonicBlockKrylovSchur(LinearOperatorBase<Field>& _Linop, GridBase* _Grid,
+                               RealD _Tolerance, ComplexD _shift = 0.0,
+                               RitzFilter _rf = EvalNormSmall)
+    : Base(_Linop, _Grid, _Tolerance, _shift, _rf)
+  {}
+
+  //--------------------------------------------------------------------
+  // Main entry point (same signature and outer structure as the base;
+  // the Schur-sort restart is replaced by harmonicRestart()).
+  //--------------------------------------------------------------------
+  virtual void operator()(const std::vector<Field>& v0, int _maxIter, int _Nm, int _Nk,
+                          int _Nstop, int _Nblock = 1, bool doubleOrthog = true,
+                          bool doVerify = false) override
+  {
+    MaxIter = _maxIter;
+    Nm      = _Nm;
+    Nk      = _Nk;
+    Nstop   = _Nstop;
+    Nblock  = _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;
+
+    RealD approxLambdaMax = approxMaxEval(v0[0]);
+    rtol = Tolerance * approxLambdaMax;
+    std::cout << GridLogMessage
+              << "TrueHarmonicBlockKrylovSchur: approx max eval = " << approxLambdaMax
+              << ", rtol = " << rtol
+              << ", shift = " << shift << std::endl;
+
+    H = CMat::Zero(N, N);
+    B = CMat::Zero(N, Nblock);
+
+    int start = 0;
+    std::vector<Field> startBlock = expandStartBlock(v0);
+
+    for (int iter = 0; iter < MaxIter; iter++) {
+      std::cout << GridLogMessage
+                << "TrueHarmonicBlockKrylovSchur: restart iteration " << iter << std::endl;
+
+      // ---- Block Arnoldi: extend from block 'start' to block Nm/Nblock ----
+      blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog);
+      std::cout << GridLogMessage << "blockArnoldiIteration done " << std::endl;
+      start = Nk/Nblock;
+
+      if (doVerify) {
+        std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi";
+        verify(lbl);
+      }
+
+      // ---- Harmonic extraction + exact thick restart ----
+      harmonicRestart();
+      std::cout << GridLogMessage << "harmonicRestart done " << std::endl;
+
+      // Restart from the residual block W (exact by identity (3))
+      startBlock = F;
+
+      if (doVerify) {
+        std::string lbl = "iter " + std::to_string(iter) + " after harmonic restart";
+        verify(lbl);
+      }
+
+      // ---- Eigensystem of retained H_k = Y^dag A Y for convergence ----
+      CMat Hk = H(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk));
+      computeEigensystem(Hk, Nk);
+
+      int Nconv = converged(Nk);
+      std::cout << GridLogMessage
+                << "TrueHarmonicBlockKrylovSchur: converged " << Nconv
+                << " / " << Nstop << std::endl;
+
+      if (Nconv >= Nstop || iter == MaxIter - 1) {
+        std::cout << GridLogMessage
+                  << "TrueHarmonicBlockKrylovSchur: done after " << iter
+                  << " restarts, " << Nconv << " converged." << std::endl;
+        std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl;
+
+        if (doEvalCheck) {
+          Field w(Grid_);
+          for (int k = 0; k < (int)evecs.size(); k++) {
+            Linop.Op(evecs[k], w);
+            ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+            w -= eval_est * evecs[k];
+            RealD res = std::sqrt(norm2(w));
+            std::cout << GridLogMessage << "TrueHarmonicBlockKrylovSchur: evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+
+        return;
+      }
+    }
+  }
+
+protected:
+
+  //--------------------------------------------------------------------
+  // Harmonic thick restart
+  //--------------------------------------------------------------------
+  // On entry:  A V = V H + F B^dag exact, V = basis (Nm fields), F Nblock
+  //            residual fields, H/B the Nm-sized small matrices.
+  // On exit:   basis = Y (Nk orthonormal fields spanning the harmonic
+  //            Ritz space), F = What (Nblock orthonormal fields, perp Y),
+  //            H top-left Nk x Nk = Y^dag A Y (dense), B top Nk rows the
+  //            new coupling; relation A Y = Y H_k + F B^dag again exact.
+  void harmonicRestart()
+  {
+    int N = Nm;
+
+    // ---- Small harmonic eigenproblem:  Hhat g = (theta - sigma) g ----
+    CMat Hs = H - shift * CMat::Identity(N, N);
+    Eigen::FullPivLU<CMat> lu(Hs.adjoint());
+    assert(lu.isInvertible() &&
+           "TrueHarmonicBlockKrylovSchur: H - shift*I singular; perturb the shift");
+    CMat M    = lu.solve(B);              // N x Nblock :  Hs^{-H} B
+    CMat Hhat = Hs + M * B.adjoint();     // N x N
+
+    Eigen::ComplexEigenSolver<CMat> es(Hhat);
+
+    // Sort harmonic values (theta - sigma) by the Ritz filter
+    // (EvalNormSmall -> closest to the shift).
+    ComplexComparator cComp(ritzFilter);
+    std::vector<int> idx(N);
+    std::iota(idx.begin(), idx.end(), 0);
+    std::sort(idx.begin(), idx.end(), [&](int a, int b){
+      return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b)));
+    });
+
+    CMat G(N, Nk);
+    CVec theta(Nk);
+    for (int k = 0; k < Nk; k++) {
+      G.col(k) = es.eigenvectors().col(idx[k]);
+      theta(k) = es.eigenvalues()(idx[k]) + toStdCmplx(shift);
+    }
+
+    std::cout << GridLogMessage
+              << "TrueHarmonicBlockKrylovSchur: harmonic Ritz values nearest shift:" << std::endl;
+    for (int k = 0; k < Nk; k++)
+      std::cout << GridLogMessage << "  [" << k << "] " << theta(k) << std::endl;
+
+    // ---- Orthonormalise the little vectors:  G = Gt Rg ----
+    Eigen::HouseholderQR<CMat> qr(G);
+    CMat Gt = qr.householderQ() * CMat::Identity(N, Nk);
+    CMat Rg = qr.matrixQR().topLeftCorner(Nk, Nk).template triangularView<Eigen::Upper>();
+    for (int k = 0; k < Nk; k++)
+      assert(std::abs(Rg(k, k)) > 1e-14 &&
+             "TrueHarmonicBlockKrylovSchur: harmonic vectors linearly dependent");
+    CMat RgInv = Rg.inverse();
+
+    // A (V G) = (V G) diag(theta) + W (B^dag G),  W = F - V M
+    //   =>  A (V Gt) = (V Gt) T + W Ct
+    CMat T  = Rg * theta.asDiagonal() * RgInv;   // Nk x Nk
+    CMat Ct = (B.adjoint() * G) * RgInv;         // Nblock x Nk
+
+    // ---- Large-field work ----
+    // Y = V Gt : new orthonormal basis (Nk fields)
+    std::vector<Field> Y;
+    constructThin(Y, basis, Gt, Nk);
+
+    // W = F - V M : Nblock residual directions containing every harmonic
+    // residual (identity (3) in the class comment)
+    std::vector<Field> W = F;
+    for (int t = 0; t < Nblock; t++)
+      for (int j = 0; j < N; j++)
+        W[t] -= basis[j] * M(j, t);
+
+    // Orthogonalise W against Y (two passes), recording P = Y^dag W
+    CMat P = CMat::Zero(Nk, Nblock);
+    for (int pass = 0; pass < 2; pass++) {
+      for (int j = 0; j < Nk; j++) {
+        for (int t = 0; t < Nblock; t++) {
+          ComplexD c = innerProduct(Y[j], W[t]);
+          P(j, t) += toStdCmplx(c);
+          W[t] -= c * Y[j];
+        }
+      }
+    }
+
+    // Block QR of the projected residual: W <- What, returns Rw
+    CMat Rw = blockQR(W);
+
+    // ---- Reassemble the exact Krylov-Schur relation ----
+    //   A Y = Y (T + P Ct) + What (Rw Ct)
+    CMat Hk = T + P * Ct;                        // Nk x Nk, dense
+    CMat Bc = Rw * Ct;                           // Nblock x Nk
+
+    H = CMat::Zero(N, N);
+    H(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)) = Hk;
+    B = CMat::Zero(N, Nblock);
+    for (int j = 0; j < Nk; j++)
+      for (int t = 0; t < Nblock; t++)
+        B(j, t) = std::conj(Bc(t, j));
+
+    basis = Y;
+    F     = W;
+    beta_k = Bc.norm();
+
+    std::cout << GridLogMessage
+              << "TrueHarmonicBlockKrylovSchur: beta_k = " << beta_k << std::endl;
+  }
+
+  //--------------------------------------------------------------------
+  // Thin basis combination: out[i] = sum_j U[j] * R(j,i), i < ncol
+  //--------------------------------------------------------------------
+  void constructThin(std::vector<Field>& out, const std::vector<Field>& U,
+                     const CMat& R, int ncol)
+  {
+    out.clear();
+    Field tmp(Grid_);
+    for (int i = 0; i < ncol; i++) {
+      tmp = Zero();
+      for (int j = 0; j < (int)U.size(); j++)
+        tmp += U[j] * R(j, i);
+      out.push_back(tmp);
+    }
+  }
+
+};
+
+NAMESPACE_END(Grid);
+
+#endif  // GRID_TRUE_HARMONIC_BLOCKED_KRYLOV_SCHUR_H

examples/Example_gamma5_block_lanczos.cc

@@ -0,0 +1,188 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./examples/Example_gamma5_block_lanczos.cc
+
+    Copyright (C) 2026
+
+    Author: Chulwoo Jung <chulwoo@bnl.gov>
+
+    γ5-Block Lanczos example for the Wilson Dirac operator.
+
+    Reads a gauge configuration from "config" (NERSC format) and Lanczos
+    parameters from "LanParams.xml".  Runs Gamma5BlockLanczos to
+    compute eigenvalues of D_W directly (not H_W = γ5 D_W).
+
+    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 <cstdlib>
+
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+
+using namespace std;
+using namespace Grid;
+
+namespace Grid {
+
+struct LanczosParameters : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
+    RealD,   mass,
+    Integer, Nstop,
+    Integer, Nk,
+    Integer, Np,
+    Integer, maxIter,
+    Integer, reorthog,
+    Integer, verify,
+    Integer, ReadEvec,
+    RealD,   resid)
+
+  LanczosParameters()
+    : mass(-0.5), Nstop(10), Nk(20), maxIter(100),
+      reorthog(1), verify(0), ReadEvec(0), resid(1e-8)
+  {}
+
+  template<class ReaderClass>
+  LanczosParameters(Reader<ReaderClass>& r) { initialize(r); }
+
+  template<class ReaderClass>
+  void initialize(Reader<ReaderClass>& r) {
+    read(r, "LanczosParameters", *this);
+  }
+};
+
+} // namespace Grid
+
+template<class T>
+void writeField(T& in, std::string const& fname) {
+#ifdef HAVE_LIME
+  std::cout << GridLogMessage << "Writing to: " << fname << std::endl;
+  Grid::emptyUserRecord record;
+  Grid::ScidacWriter WR(in.Grid()->IsBoss());
+  WR.open(fname);
+  WR.writeScidacFieldRecord(in, record, 0);
+  WR.close();
+#endif
+}
+
+typedef WilsonFermionD WilsonOp;
+typedef typename WilsonFermionD::FermionField FermionField;
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
+      GridDefaultLatt(),
+      GridDefaultSimd(Nd, vComplex::Nsimd()),
+      GridDefaultMpi());
+  GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  std::vector<int> seeds4({1, 2, 3, 4});
+  GridParallelRNG RNG4(UGrid);
+  RNG4.SeedFixedIntegers(seeds4);
+
+  // Read gauge configuration
+  LatticeGaugeField Umu(UGrid);
+  FieldMetaData header;
+  std::string configFile("config");
+  NerscIO::readConfiguration(Umu, header, configFile);
+  std::cout << GridLogMessage << "Loaded gauge configuration: " << configFile << std::endl;
+
+  // Read Lanczos parameters
+  LanczosParameters Params;
+  {
+    XmlReader rd("LanParams.xml");
+    read(rd, "LanczosParameters", Params);
+  }
+  std::cout << GridLogMessage << Params << std::endl;
+
+  // Build Wilson Dirac operator and wrap in a non-Hermitian linear operator
+  std::vector<Complex> boundary = {1, 1, 1, -1};
+  WilsonOp::ImplParams WilsonParams(boundary);
+  WilsonOp Dwilson(Umu, *UGrid, *UrbGrid, Params.mass, WilsonParams);
+  NonHermitianLinearOperator<WilsonOp, FermionField> DLinOp(Dwilson);
+
+  // γ5 functor: for 4D Wilson fermions γ5 is Gamma(Gamma5)
+  Gamma G5(Gamma::Algebra::Gamma5);
+  auto gamma5 = [&G5](const FermionField& in, FermionField& out) {
+    out = G5 * in;
+  };
+
+  // Starting vectors: two independent random vectors
+  FermionField src(UGrid), src2(UGrid);
+  random(RNG4, src);
+  random(RNG4, src2);
+  std::cout << GridLogMessage << "Using two random starting vectors" << std::endl;
+
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "*******************************************" << std::endl;
+  std::cout << GridLogMessage << " Running γ5-Block Lanczos" << std::endl;
+  std::cout << GridLogMessage << "   mass     = " << Params.mass     << std::endl;
+  std::cout << GridLogMessage << "   Nk       = " << Params.Nk       << std::endl;
+  std::cout << GridLogMessage << "   maxIter  = " << Params.maxIter  << std::endl;
+  std::cout << GridLogMessage << "   Nstop    = " << Params.Nstop    << std::endl;
+  std::cout << GridLogMessage << "   resid    = " << Params.resid    << std::endl;
+  std::cout << GridLogMessage << "   reorthog = " << Params.reorthog << std::endl;
+  std::cout << GridLogMessage << "   verify   = " << Params.verify   << std::endl;
+  std::cout << GridLogMessage << "*******************************************" << std::endl;
+  std::cout << GridLogMessage << std::endl;
+
+  Gamma5BlockLanczos<FermionField> G5BL(DLinOp, UGrid, gamma5, Params.resid);
+  G5BL.doEvalCheck = (Params.verify != 0);
+  G5BL.doVerify = (Params.verify != 0);
+//  G5BL(src, src2, Params.maxIter, Params.Nstop, Params.reorthog != 0);
+  G5BL.restart(src, src2, Params.maxIter, Params.Nk+Params.Np, Params.Nk, Params.Nstop, Params.reorthog != 0, EvalNormSmall);
+//  G5BL.implicitRestart(src, src2, Params.maxIter, Params.Nk+Params.Np, Params.Nk, Params.Nstop, Params.reorthog != 0, EvalNormSmall);
+  if (Params.verify ) G5BL.verify("after restart");
+
+  // Summary of results
+  Eigen::VectorXcd        evals     = G5BL.getEvals();
+  std::vector<RealD>      residuals = G5BL.getResiduals();
+  std::vector<FermionField> evecs   = G5BL.getEvecs();
+
+  int Nout = (int)evals.size();
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "*******************************************" << std::endl;
+  std::cout << GridLogMessage << " γ5-Block Lanczos: " << Nout << " Ritz pairs" << std::endl;
+  std::cout << GridLogMessage << "*******************************************" << std::endl;
+  for (int i = 0; i < Nout; i++) {
+    std::cout << GridLogMessage
+              << "  [" << std::setw(3) << i << "]"
+              << "  lambda = " << evals(i)
+              << "  |res| = " << residuals[i] << std::endl;
+  }
+
+  // Write the first Nstop eigenvectors (in SCIDAC format when LIME is available)
+  int Nwrite = std::min((int)Params.Nstop, Nout);
+  for (int i = 0; i < Nwrite; i++) {
+    std::string fname = "./g5bl_evec_m" + std::to_string(Params.mass)
+                        + "_" + std::to_string(i);
+    writeField(evecs[i], fname);
+  }
+
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "Done" << std::endl;
+
+  Grid_finalize();
+  return 0;
+}

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
@@ -52,6 +54,9 @@ struct LanczosParameters: Serializable {
                                 Integer, Np,
                                 Integer, ReadEvec,
                                 Integer, maxIter,
+                                Integer, Nblock,
+                                Integer, verify,
+		  		RealD, shift , 
 	  			RealD, resid,
 	  			RealD, ChebyLow,
 	  			RealD, ChebyHigh,
@@ -113,12 +118,9 @@ public:
   InvertNonHermitianLinearOperator(Matrix &Mat,RealD stp=1e-8): _Mat(Mat),_stp(stp){};
   // Support for coarsening to a multigrid
   void OpDiag (const Field &in, Field &out) {
-//    _Mat.Mdiag(in,out);
-//    out = out + shift*in;
     assert(0);
   }
   void OpDir  (const Field &in, Field &out,int dir,int disp) {
-//    _Mat.Mdir(in,out,dir,disp);
     assert(0);
   }
   void OpDirAll  (const Field &in, std::vector<Field> &out){
@@ -127,11 +129,6 @@ public:
   };
   void Op     (const Field &in, Field &out){
     Field tmp(in.Grid());
-//    _Mat.M(in,out);
-//  RealD mass=-shift;
-//  WilsonCloverFermionD Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t);
-//  NonHermitianLinearOperator<Matrix,Field> HermOp(_Mat);
-//  BiCGSTAB<Field> CG(_stp,10000);
     _Mat.Mdag(in,tmp);
     MdagMLinearOperator<Matrix,Field> HermOp(_Mat);
     ConjugateGradient<Field> CG(_stp,10000);
@@ -309,12 +306,21 @@ 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";
     std::cout << GridLogIRL<< "Reading evecs from "<<evecs_file<<std::endl;
-    emptyUserRecord record;
+    Grid::emptyUserRecord record;
     Grid::ScidacReader RD;
     RD.open(evecs_file);
     RD.readScidacFieldRecord(src[0],record);
@@ -326,29 +332,68 @@ int main (int argc, char ** argv)
   std::cout << "[DEBUG] Source at origin = " <<  tmpSrc << std::endl;
   LatticeFermion src2 = src[0];
 
-  // Run KrylovSchur and Arnoldi on a Hermitian matrix
-  std::cout << GridLogMessage << "Running Krylov Schur" << std::endl;
-  // KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge);
-//  KrylovSchur KrySchur (Dsq, FGrid, 1e-8,EvalImNormSmall);
-//  KrySchur(src, maxIter, Nm, Nk, Nstop);
-//  KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall);
-//  Hacked, really EvalImagSmall
-#if 1
-    RealD shift=1.5;
-    KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
-    KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
-#else
-    KrylovSchur KrySchur (Iwilson, UGrid, resid,EvalImNormSmall);
-    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
-#endif
-  std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl;
+  // Run KrylovSchur / BlockKrylovSchur / TrueHarmonicBlockKrylovSchur on a Hermitian matrix
+    RealD shift=LanParams.shift;
 
-  src[0]=KrySchur.evecs[0];
-  for (int i=1;i<Nstop;i++) src[0]+=KrySchur.evecs[i];
-  for (int i=0;i<Nstop;i++) 
+    // --use-blockKS : off -> plain (non-block) KrylovSchur
+    //                 on  -> block family (BlockKrylovSchur / TrueHarmonicBlockKrylovSchur)
+    // --use-harmonic: off -> no shift/harmonic Ritz extraction
+    //                 on, without --use-blockKS -> plain KrylovSchur with shifted (harmonic) Schur extraction
+    //                 on, with --use-blockKS    -> TrueHarmonicBlockKrylovSchur (true harmonic Ritz)
+    bool useBlockKS  = GridCmdOptionExists(argv, argv+argc, std::string("--use-blockKS"));
+    bool useHarmonic = GridCmdOptionExists(argv, argv+argc, std::string("--use-harmonic"));
+
+    std::vector<LatticeFermion> finalEvecs;
+
+    if (!useBlockKS) {
+      KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+      KrySchur.doEvalCheck=true;
+      if (useHarmonic) {
+        std::cout << GridLogMessage << "Running KrylovSchur (shifted/harmonic Ritz)" << std::endl;
+        KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
+      } else {
+        std::cout << GridLogMessage << "Running KrylovSchur" << std::endl;
+        KrySchur(src[0], maxIter, Nm, Nk, Nstop);
+      }
+      std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
+      finalEvecs = KrySchur.evecs;
+    } else {
+      int Nblock=4;
+      Nblock=LanParams.Nblock;
+      bool if_verify=false;
+      if(LanParams.verify) if_verify=true;
+
+      auto gamma5Lambda = [](const FermionField& v, FermionField& out) {
+        Gamma g5(Gamma::Algebra::Gamma5);
+        out = g5 * v;
+      };
+
+      if (useHarmonic) {
+        std::cout << GridLogMessage << "Running TrueHarmonicBlockKrylovSchur" << std::endl;
+        TrueHarmonicBlockKrylovSchur<FermionField> KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall);
+        KrySchur.doEvalCheck=true;
+        KrySchur.gamma5Func = gamma5Lambda;
+        KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify);
+        std::cout << GridLogMessage << "TrueHarmonicBlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
+        finalEvecs = KrySchur.evecs;
+      } else {
+        std::cout << GridLogMessage << "Running BlockKrylovSchur" << std::endl;
+        BlockKrylovSchur<FermionField> KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+        KrySchur.doEvalCheck=true;
+        KrySchur.gamma5Func = gamma5Lambda;
+        KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify);
+        std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
+        finalEvecs = KrySchur.evecs;
+      }
+    }
+
+  src[0]=finalEvecs[0];
+  std::cout << GridLogMessage << "finalEvecs= "<< finalEvecs.size() <<std::endl;
+  for (int i=1;i<Nk;i++) src[0]+=finalEvecs[i];
+  for (int i=0;i<Nstop;i++)
   {
 	std::string evfile ("./evec_"+std::to_string(mass)+"_"+std::to_string(i));
-        auto evdensity = localInnerProduct(KrySchur.evecs[i],KrySchur.evecs[i] );
+        auto evdensity = localInnerProduct(finalEvecs[i],finalEvecs[i] );
         writeFile(evdensity,evfile);
 
   }

examples/LanParams.xml

@@ -0,0 +1,22 @@
+<?xml version="1.0"?>
+<grid>
+  <LanczosParameters>
+    <mass>0</mass>
+    <mstep>-0.025</mstep>
+    <M5>1.8</M5>
+    <Ls>48</Ls>
+    <Nstop>100</Nstop>
+    <Nk>100</Nk>
+    <Np>100</Np>
+    <ReadEvec>0</ReadEvec>
+    <maxIter>1000</maxIter>
+    <reorthog>1</reorthog>
+    <Nblock>4</Nblock>
+    <verify>0</verify>
+    <shift>1.5</shift>
+    <resid>1e-8</resid>
+    <ChebyLow>1</ChebyLow>
+    <ChebyHigh>100</ChebyHigh>
+    <ChebyOrder>51</ChebyOrder>
+  </LanczosParameters>
+</grid>

examples/fft5d.cc

@@ -0,0 +1,469 @@
+/*************************************************************************************
+
+  fft5d.cc — Fourier analysis of a time series of 4-D lattice scalar fields.
+
+  Designed for force-norm files from FTHMC (one RealD per site, SCIDAC or binary).
+  Assembles a (4+1)-D structure [trajectory][t][z][y][x] using Grid lattice fields
+  and performs three FFT analyses:
+
+    --fft spatial
+        4-D spatial FFT using Grid's FFT class (MPI+GPU parallel via Cshift/cufft).
+        Outputs shell-averaged P(|k|^2) averaged over trajectories, and a
+        per-mode table P(kt,kz,ky,kx).
+
+    --fft traj
+        1-D trajectory-axis FFT at each lattice site using FFTW3 locally at each
+        MPI rank.  Outputs site-averaged P(f_traj) and autocorrelation C(lag).
+        Results are combined with GlobalSumVector.
+
+    --fft all
+        Both of the above, plus a 2-D cross spectrum P(f_traj, |k_spatial|^2):
+        spatial FFT per trajectory, then trajectory FFTW on the momentum-space series.
+
+  Normalisations:
+    Spatial FFT  (volume V): P_spatial(k) = |F(k)|^2 / V^2
+      => (1/V) * sum_k P(k) = site mean-square per trajectory
+    Trajectory FFT (Ntraj):  P_traj(f) = |F(f)|^2 / Ntraj^2
+      => sum_f P(f) = site mean-square over trajectories
+
+  Build: add to examples/Make.inc (see bottom of this file), then make.
+
+  Usage (follows Grid conventions):
+    fft5d --grid Lx.Ly.Lz.Lt [--mpi Px.Py.Pz.Pt] \
+          [--fft spatial|traj|all] [--format scidac|binary] \
+          [--output PREFIX] file1 file2 ...
+
+*************************************************************************************/
+
+#include <Grid/Grid.h>
+// Grid's FFT.h uses cufft on CUDA builds; for the trajectory axis we also need
+// CPU-side FFTW3 (already linked via -lfftw3 in GRID_LIBS).
+#include <fftw3.h>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+
+using namespace Grid;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Global coordinate of a local site (osite, lane) on this MPI rank
+// ─────────────────────────────────────────────────────────────────────────────
+static void globalCoor(int osite, int lane, GridCartesian* g, Coordinate& gc)
+{
+    Coordinate oc, ic;
+    Lexicographic::CoorFromIndex(oc, osite, g->_rdimensions);
+    Lexicographic::CoorFromIndex(ic, lane,  g->_simd_layout);
+    for (int d = 0; d < Nd; d++)
+        gc[d] = g->_processor_coor[d] * g->_ldimensions[d]
+              + oc[d] * g->_simd_layout[d] + ic[d];
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Shell-average a power LatticeRealD over |k|^2 bins (MPI-aware via GlobalSumVector)
+// ─────────────────────────────────────────────────────────────────────────────
+static void shellAverage(const LatticeRealD& power, GridCartesian* g,
+                         double norm, std::ostream& ofs)
+{
+    const int Nsimd = vRealD::Nsimd();
+    Coordinate fdims = g->_fdimensions;
+    int maxk2 = 0;
+    for (int d = 0; d < Nd; d++) { int h = fdims[d]/2; maxk2 += h*h; }
+
+    std::vector<RealD> psum(maxk2+1, 0.0), cnt(maxk2+1, 0.0);
+    {
+        auto pv = power.View(CpuRead);
+        Coordinate gc(Nd);
+        for (int os = 0; os < (int)g->oSites(); os++) {
+            for (int lane = 0; lane < Nsimd; lane++) {
+                globalCoor(os, lane, g, gc);
+                int k2 = 0;
+                for (int d = 0; d < Nd; d++) {
+                    int kd = std::min(gc[d], fdims[d] - gc[d]);
+                    k2 += kd*kd;
+                }
+                psum[k2] += (RealD)extractLane(lane, pv[os]);
+                cnt [k2] += 1.0;
+            }
+        }
+    }
+    g->GlobalSumVector(psum.data(), (int)psum.size());
+    g->GlobalSumVector(cnt .data(), (int)cnt .size());
+
+    for (int k2 = 0; k2 <= maxk2; k2++)
+        if (cnt[k2] > 0)
+            ofs << k2 << "  " << std::sqrt((double)k2) << "  "
+                << (psum[k2] / cnt[k2]) * norm << "\n";
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Convert LatticeRealD → LatticeComplexD (zero imaginary part)
+// ─────────────────────────────────────────────────────────────────────────────
+static LatticeComplexD toComplex(const LatticeRealD& r)
+{
+    std::vector<RealD>    lr;
+    unvectorizeToLexOrdArray(lr, r);
+    std::vector<ComplexD> lc(lr.size());
+    for (size_t i = 0; i < lr.size(); i++) lc[i] = ComplexD(lr[i], 0.0);
+    LatticeComplexD c(r.Grid());
+    vectorizeFromLexOrdArray(lc, c);
+    return c;
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// File readers
+// ─────────────────────────────────────────────────────────────────────────────
+static LatticeComplexD readScidac(const std::string& fname, GridCartesian* g)
+{
+    LatticeRealD field(g);
+    emptyUserRecord rec;
+    ScidacReader RD;
+    RD.open(fname);
+    RD.readScidacFieldRecord(field, rec);
+    RD.close();
+    return toComplex(field);
+}
+
+static LatticeComplexD readBinary(const std::string& fname, GridCartesian* g)
+{
+    // Raw IEEE doubles in lex order [x][y][z][t] written serially.
+    // Rank 0 reads the file, broadcasts to all ranks via GlobalSumVector.
+    Coordinate fdims = g->_fdimensions;
+    long vol4 = 1; for (int d = 0; d < Nd; d++) vol4 *= fdims[d];
+
+    std::vector<RealD> buf(vol4, 0.0);
+    if (g->IsBoss()) {
+        std::ifstream f(fname, std::ios::binary);
+        if (!f) throw std::runtime_error("Cannot open: " + fname);
+        f.read(reinterpret_cast<char*>(buf.data()), vol4 * sizeof(RealD));
+    }
+    g->GlobalSumVector(buf.data(), (int)vol4);   // broadcast from rank 0
+
+    LatticeRealD field(g);
+    Coordinate ldims = g->_ldimensions, pcoor = g->_processor_coor;
+    for (long ls = 0; ls < g->lSites(); ls++) {
+        Coordinate lc;
+        g->LocalIndexToLocalCoor(ls, lc);
+        long glex = 0, stride = 1;
+        for (int d = 0; d < Nd; d++) {
+            glex += (pcoor[d]*ldims[d] + lc[d]) * stride;
+            stride *= fdims[d];
+        }
+        pokeLocalSite(buf[glex], field, lc);
+    }
+    return toComplex(field);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Load all files into a trajectory vector
+// ─────────────────────────────────────────────────────────────────────────────
+static std::vector<LatticeComplexD>
+loadFiles(const std::vector<std::string>& files, GridCartesian* g,
+          const std::string& fmt)
+{
+    std::vector<LatticeComplexD> traj;
+    traj.reserve(files.size());
+    for (int n = 0; n < (int)files.size(); n++) {
+        std::cout << GridLogMessage << "[" << n << "] reading " << files[n] << "\n";
+        if (fmt == "scidac") traj.push_back(readScidac(files[n], g));
+        else                 traj.push_back(readBinary (files[n], g));
+    }
+    return traj;
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Spatial FFT analysis
+// ─────────────────────────────────────────────────────────────────────────────
+static void analyzeSpatial(const std::vector<LatticeComplexD>& traj,
+                            GridCartesian* g, const std::string& pfx)
+{
+    int  Ntraj = (int)traj.size();
+    long vol4  = 1; for (int d = 0; d < Nd; d++) vol4 *= g->_fdimensions[d];
+
+    FFT theFFT(g);
+    LatticeRealD pavg(g); pavg = Zero();
+    for (int n = 0; n < Ntraj; n++) {
+        LatticeComplexD fk(g);
+        theFFT.FFT_all_dim(fk, traj[n], FFT::forward);
+        // Evaluate product before applying real() — real() is not defined for
+        // unevaluated LatticeBinaryExpression.
+        LatticeComplexD fk_sq(g); fk_sq = conjugate(fk) * fk;
+        std::vector<ComplexD> lc; unvectorizeToLexOrdArray(lc, fk_sq);
+        std::vector<RealD> lr(lc.size());
+        for (size_t i = 0; i < lc.size(); i++) lr[i] = lc[i].real();
+        LatticeRealD pk(g); vectorizeFromLexOrdArray(lr, pk);
+        pavg += pk;
+    }
+    pavg *= (1.0 / Ntraj);
+
+    // Shell-averaged spectrum
+    {
+        std::ofstream fs(pfx + "_spatial_shell.dat");
+        fs << "# k2  |k|  P_shell_avg   (P = |F|^2 / V^2 / shell_count)\n";
+        shellAverage(pavg, g, 1.0 / ((double)vol4 * vol4), fs);
+        std::cout << GridLogMessage << "Written: " << pfx << "_spatial_shell.dat\n";
+    }
+
+    // Per-mode table (rank 0 only, via peekSite)
+    if (g->IsBoss()) {
+        std::ofstream fm(pfx + "_spatial_modes.dat");
+        fm << "# kt  kz  ky  kx  k2  |k|  P_traj_avg\n";
+        Coordinate fdims = g->_fdimensions, site(Nd);
+        double norm = 1.0 / ((double)vol4 * vol4);
+        for (site[3]=0; site[3]<fdims[3]; site[3]++) { int kt=std::min(site[3],fdims[3]-site[3]);
+        for (site[2]=0; site[2]<fdims[2]; site[2]++) { int kz=std::min(site[2],fdims[2]-site[2]);
+        for (site[1]=0; site[1]<fdims[1]; site[1]++) { int ky=std::min(site[1],fdims[1]-site[1]);
+        for (site[0]=0; site[0]<fdims[0]; site[0]++) { int kx=std::min(site[0],fdims[0]-site[0]);
+            RealD p; peekSite(p, pavg, site);
+            int k2 = kt*kt+kz*kz+ky*ky+kx*kx;
+            fm << kt<<" "<<kz<<" "<<ky<<" "<<kx<<" "<<k2<<" "
+               <<std::sqrt((double)k2)<<" "<<p*norm<<"\n";
+        }}}}
+        std::cout << GridLogMessage << "Written: " << pfx << "_spatial_modes.dat\n";
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Local FFTW trajectory FFT: each rank FFTs its own local sites independently
+// ─────────────────────────────────────────────────────────────────────────────
+static void trajFFT(const std::vector<LatticeComplexD>& in,
+                    std::vector<LatticeComplexD>& out)
+{
+    int    Ntraj  = (int)in.size();
+    GridBase* g   = in[0].Grid();
+    long   lsites = g->lSites();
+
+    // fftw_complex is double[2] — std::vector<double[2]> is not supported by nvcc.
+    // Use std::vector<double> with 2x size and reinterpret_cast.
+    std::vector<double> ibuf((long)Ntraj * lsites * 2, 0.0);
+    for (int n = 0; n < Ntraj; n++) {
+        std::vector<ComplexD> lc;
+        unvectorizeToLexOrdArray(lc, in[n]);
+        for (long s = 0; s < lsites; s++) {
+            ibuf[((long)n*lsites + s)*2 + 0] = lc[s].real();
+            ibuf[((long)n*lsites + s)*2 + 1] = lc[s].imag();
+        }
+    }
+
+    // lsites transforms of length Ntraj, stride=lsites, dist=1
+    std::vector<double> obuf((long)Ntraj * lsites * 2, 0.0);
+    fftw_complex* iptr = reinterpret_cast<fftw_complex*>(ibuf.data());
+    fftw_complex* optr = reinterpret_cast<fftw_complex*>(obuf.data());
+    int n1[1] = {Ntraj};
+    fftw_plan p = fftw_plan_many_dft(
+        1, n1, (int)lsites,
+        iptr, nullptr, (int)lsites, 1,
+        optr, nullptr, (int)lsites, 1,
+        FFTW_FORWARD, FFTW_ESTIMATE);
+    fftw_execute(p);
+    fftw_destroy_plan(p);
+
+    // vector::assign triggers _M_fill_assign which needs a default constructor;
+    // Lattice has none.  Use explicit emplace_back instead.
+    out.clear(); out.reserve(Ntraj);
+    for (int k = 0; k < Ntraj; k++) out.emplace_back(g);
+    for (int k = 0; k < Ntraj; k++) {
+        std::vector<ComplexD> lc(lsites);
+        for (long s = 0; s < lsites; s++)
+            lc[s] = ComplexD(obuf[((long)k*lsites+s)*2], obuf[((long)k*lsites+s)*2+1]);
+        vectorizeFromLexOrdArray(lc, out[k]);
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Trajectory-axis FFT analysis
+// ─────────────────────────────────────────────────────────────────────────────
+static void analyzeTraj(const std::vector<LatticeComplexD>& traj,
+                         GridCartesian* g, const std::string& pfx)
+{
+    int  Ntraj = (int)traj.size();
+    int  Nf    = Ntraj / 2 + 1;
+    long vol4  = 1; for (int d = 0; d < Nd; d++) vol4 *= g->_fdimensions[d];
+
+    std::vector<LatticeComplexD> ftraj;
+    trajFFT(traj, ftraj);
+
+    // P(k) = (1/vol4) * sum_sites |F_traj(k)|^2 / Ntraj^2
+    std::vector<RealD> Pavg(Ntraj, 0.0);
+    for (int k = 0; k < Ntraj; k++) {
+        LatticeComplexD tmp(g); tmp = conjugate(ftraj[k]) * ftraj[k];
+        std::vector<ComplexD> lc; unvectorizeToLexOrdArray(lc, tmp);
+        std::vector<RealD> lr(lc.size());
+        for (size_t i = 0; i < lc.size(); i++) lr[i] = lc[i].real();
+        LatticeRealD pk(g); vectorizeFromLexOrdArray(lr, pk);
+        RealD s = 0.0;
+        for (long ls = 0; ls < g->lSites(); ls++) {
+            Coordinate lc; g->LocalIndexToLocalCoor(ls, lc);
+            RealD v; peekLocalSite(v, pk, lc);
+            s += v;
+        }
+        g->GlobalSum(s);
+        Pavg[k] = s / ((double)vol4 * Ntraj * Ntraj);
+    }
+
+    {
+        std::ofstream f(pfx + "_traj_power.dat");
+        f << "# k  freq  P_avg   (sum_k P = site mean-square)\n";
+        for (int k = 0; k < Nf; k++)
+            f << k << "  " << (double)k/Ntraj << "  " << Pavg[k] << "\n";
+        std::cout << GridLogMessage << "Written: " << pfx << "_traj_power.dat\n";
+    }
+
+    // Autocorrelation via IFFT of the power spectrum
+    std::vector<double> Pc_buf(Nf * 2, 0.0);
+    for (int k = 0; k < Nf; k++) { Pc_buf[k*2] = Pavg[k]; Pc_buf[k*2+1] = 0.0; }
+    fftw_complex* Pc = reinterpret_cast<fftw_complex*>(Pc_buf.data());
+    std::vector<double> acorr(Ntraj, 0.0);
+    fftw_plan ip = fftw_plan_dft_c2r(1, &Ntraj, Pc, acorr.data(), FFTW_ESTIMATE);
+    fftw_execute(ip); fftw_destroy_plan(ip);
+    double c0 = acorr[0] / Ntraj;
+    {
+        std::ofstream f(pfx + "_traj_autocorr.dat");
+        f << "# lag  C(lag)  C(lag)/C(0)\n";
+        for (int lag = 0; lag < Ntraj/2; lag++) {
+            double c = acorr[lag] / Ntraj;
+            f << lag << "  " << c << "  " << (c0 > 0 ? c/c0 : 0.0) << "\n";
+        }
+        std::cout << GridLogMessage << "Written: " << pfx << "_traj_autocorr.dat\n";
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Full 5-D analysis: spatial FFT → trajectory FFT → 2-D cross spectrum
+// ─────────────────────────────────────────────────────────────────────────────
+static void analyzeAll5D(const std::vector<LatticeComplexD>& traj,
+                          GridCartesian* g, const std::string& pfx)
+{
+    int  Ntraj = (int)traj.size();
+    int  Nf    = Ntraj / 2 + 1;
+    long vol4  = 1; for (int d = 0; d < Nd; d++) vol4 *= g->_fdimensions[d];
+    Coordinate fdims = g->_fdimensions;
+
+    // Step 1: spatial FFT for each trajectory
+    FFT theFFT(g);
+    std::vector<LatticeComplexD> sfft(Ntraj, LatticeComplexD(g));
+    for (int n = 0; n < Ntraj; n++)
+        theFFT.FFT_all_dim(sfft[n], traj[n], FFT::forward);
+
+    // Step 2: trajectory FFTW on the momentum-space series
+    std::vector<LatticeComplexD> tfft;
+    trajFFT(sfft, tfft);
+
+    // Step 3: P(f_traj, |k_spatial|^2) shell-averaged
+    int maxk2 = 0;
+    for (int d = 0; d < Nd; d++) { int h = fdims[d]/2; maxk2 += h*h; }
+
+    long nb = (long)Nf * (maxk2+1);
+    std::vector<RealD> p2d(nb, 0.0), cnt2d(nb, 0.0);
+    const int Nsimd = vComplexD::Nsimd();
+    double norm = 1.0 / ((double)Ntraj * Ntraj * vol4 * vol4);
+
+    for (int kf = 0; kf < Nf; kf++) {
+        LatticeComplexD tmp(g); tmp = conjugate(tfft[kf]) * tfft[kf];
+        std::vector<ComplexD> lc2; unvectorizeToLexOrdArray(lc2, tmp);
+        std::vector<RealD> lr2(lc2.size());
+        for (size_t i = 0; i < lc2.size(); i++) lr2[i] = lc2[i].real();
+        LatticeRealD pk(g); vectorizeFromLexOrdArray(lr2, pk);
+        auto pv = pk.View(CpuRead);
+        Coordinate gc(Nd);
+        for (int os = 0; os < (int)g->oSites(); os++) {
+            for (int lane = 0; lane < Nsimd; lane++) {
+                globalCoor(os, lane, g, gc);
+                int k2 = 0;
+                for (int d = 0; d < Nd; d++) {
+                    int kd = std::min(gc[d], fdims[d] - gc[d]);
+                    k2 += kd*kd;
+                }
+                long b = (long)kf*(maxk2+1) + k2;
+                p2d [b] += (RealD)extractLane(lane, pv[os]);
+                cnt2d[b] += 1.0;
+            }
+        }
+    }
+    g->GlobalSumVector(p2d .data(), (int)nb);
+    g->GlobalSumVector(cnt2d.data(), (int)nb);
+
+    if (g->IsBoss()) {
+        std::ofstream f(pfx + "_5d_cross.dat");
+        f << "# k_traj  freq_traj  k2_spatial  |k_spatial|  P_avg\n";
+        for (int kf = 0; kf < Nf; kf++) {
+            double freq = (double)kf / Ntraj;
+            for (int k2 = 0; k2 <= maxk2; k2++) {
+                long b = (long)kf*(maxk2+1) + k2;
+                if (cnt2d[b] > 0)
+                    f << kf << "  " << freq << "  " << k2 << "  "
+                      << std::sqrt((double)k2) << "  "
+                      << (p2d[b]/cnt2d[b]) * norm << "\n";
+            }
+            f << "\n";   // blank line between kf slices for gnuplot pm3d
+        }
+        std::cout << GridLogMessage << "Written: " << pfx << "_5d_cross.dat\n";
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Main
+// ─────────────────────────────────────────────────────────────────────────────
+static void usage(const char* prog)
+{
+    std::cerr
+        << "Usage: " << prog << " --grid Lx.Ly.Lz.Lt [--mpi Px.Py.Pz.Pt]\n"
+        << "       [--fft spatial|traj|all] [--format scidac|binary]\n"
+        << "       [--output PREFIX] file1 [file2 ...]\n";
+    exit(1);
+}
+
+int main(int argc, char** argv)
+{
+    Grid_init(&argc, &argv);   // consumes --grid, --mpi, etc.
+
+    std::string fftMode = "all";
+    std::string fmt     = "scidac";
+    std::string pfx     = "fft5d";
+    std::vector<std::string> files;
+
+    for (int i = 1; i < argc; i++) {
+        std::string a = argv[i];
+        if      (a == "--fft"    && i+1 < argc) fftMode = argv[++i];
+        else if (a == "--format" && i+1 < argc) fmt     = argv[++i];
+        else if (a == "--output" && i+1 < argc) pfx     = argv[++i];
+        else if (!a.empty() && a[0] != '-')     files.push_back(a);
+        // unknown --flags skipped (may be Grid flags already consumed)
+    }
+
+    if (files.empty()) usage(argv[0]);
+
+    GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(),
+        GridDefaultSimd(Nd, vComplexD::Nsimd()),
+        GridDefaultMpi());
+
+    Coordinate latt = GridDefaultLatt();
+    std::cout << GridLogMessage << "Lattice : "
+              << latt[0]<<"x"<<latt[1]<<"x"<<latt[2]<<"x"<<latt[3] << "\n"
+              << GridLogMessage << "Ntraj   : " << files.size() << "\n"
+              << GridLogMessage << "FFT     : " << fftMode << "\n"
+              << GridLogMessage << "Format  : " << fmt << "\n"
+              << GridLogMessage << "Output  : " << pfx << "\n";
+
+    auto traj = loadFiles(files, grid, fmt);
+
+    if (fftMode == "spatial" || fftMode == "all") analyzeSpatial(traj, grid, pfx);
+    if (fftMode == "traj"    || fftMode == "all") analyzeTraj   (traj, grid, pfx);
+    if (fftMode == "all")                         analyzeAll5D  (traj, grid, pfx);
+
+    Grid_finalize();
+    return 0;
+}
+
+/*
+  Add to examples/Make.inc to build:
+
+    bin_PROGRAMS += ... fft5d
+    fft5d_SOURCES = fft5d.cc
+    fft5d_LDADD   = $(top_builddir)/Grid/libGrid.a
+*/

tests/core/Test_plaquette_stats.cc

@@ -0,0 +1,271 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: tests/core/Test_plaquette_stats.cc
+
+Copyright (C) 2015
+
+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 */
+
+/**
+ * Test_plaquette_stats
+ *
+ * Measure every plaquette Re Tr[U_mu(x) U_nu(x+mu) U_mu†(x+nu) U_nu†(x)] / Nc
+ * across all sites and all (mu,nu) planes and report max, min, and average.
+ *
+ * Usage:
+ *   ./Test_plaquette_stats [Grid options] [--file <nersc_config>] [--hot]
+ *                          [--threshold <val>]
+ *
+ *   --file <path>       Read gauge field from NERSC-format file
+ *   --hot               Use random (hot) SU(3) start (default: cold/unit start)
+ *   --threshold <val>   Print coordinates of every plaquette with Re Tr/Nc < val
+ *
+ * Grid size defaults to 8^3 x 16; override with --grid (e.g. --grid 4.4.4.8).
+ */
+
+#include <Grid/Grid.h>
+#include <limits>
+
+using namespace std;
+using namespace Grid;
+
+template<class T>
+static void writeFile(T& in, const std::string& fname)
+{
+#ifdef HAVE_LIME
+  std::cout << GridLogMessage << "Writing: " << fname << std::endl;
+  Grid::emptyUserRecord record;
+  Grid::ScidacWriter WR(in.Grid()->IsBoss());
+  WR.open(fname);
+  WR.writeScidacFieldRecord(in, record, 0);
+  WR.close();
+#endif
+}
+
+// Return the plane label string for output
+static std::string planeName(int mu, int nu)
+{
+  const char dirs[] = "xyzt";
+  std::string s;
+  s += dirs[mu];
+  s += dirs[nu];
+  return s;
+}
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  // ---- Lattice geometry ----
+  Coordinate latt_size   = GridDefaultLatt();
+  if (latt_size.size() == 0) {
+    latt_size = Coordinate(std::vector<int>{8, 8, 8, 16});
+  }
+  Coordinate simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+
+  GridCartesian grid(latt_size, simd_layout, mpi_layout);
+
+  std::cout << GridLogMessage << "Lattice: ";
+  for (int d = 0; d < Nd; d++)
+    std::cout << latt_size[d] << (d < Nd-1 ? "x" : "\n");
+
+  // ---- Gauge field ----
+  LatticeGaugeField Umu(&grid);
+
+  // Check for --file argument
+  std::string config_file = "";
+  for (int i = 1; i < argc - 1; i++) {
+    if (std::string(argv[i]) == "--file") {
+      config_file = argv[i+1];
+      break;
+    }
+  }
+
+  bool doHot = false;
+  for (int i = 1; i < argc; i++) {
+    if (std::string(argv[i]) == "--hot") { doHot = true; break; }
+  }
+
+  double threshold = std::numeric_limits<double>::quiet_NaN();
+  bool   doThresh  = false;
+  for (int i = 1; i < argc - 1; i++) {
+    if (std::string(argv[i]) == "--threshold") {
+      threshold = std::stod(argv[i+1]);
+      doThresh  = true;
+      break;
+    }
+  }
+
+  if (!config_file.empty()) {
+    std::cout << GridLogMessage << "Reading gauge field from " << config_file << std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu, header, config_file);
+  } else {
+    std::vector<int> seeds({1, 2, 3, 4});
+    GridParallelRNG pRNG(&grid);
+    pRNG.SeedFixedIntegers(seeds);
+    if (doHot) {
+      std::cout << GridLogMessage << "Generating hot (random SU(3)) start" << std::endl;
+      SU<Nc>::HotConfiguration(pRNG, Umu);
+    } else {
+      std::cout << GridLogMessage << "Using cold (unit) gauge start" << std::endl;
+      SU<Nc>::ColdConfiguration(pRNG, Umu);
+    }
+  }
+
+  // ---- Extract link matrices ----
+  std::vector<LatticeColourMatrix> U(Nd, &grid);
+  for (int mu = 0; mu < Nd; mu++)
+    U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+
+  // ---- Per-plane plaquette statistics ----
+  //
+  // For each (mu, nu) plane (mu > nu) compute
+  //   P_munu(x) = Re Tr[plaquette] / Nc
+  // then report max, min, mean over all sites.
+  //
+  // WilsonLoops::traceDirPlaquette gives Tr[U_mu U_nu U_mu† U_nu†] (complex).
+
+  double vol = grid.gSites();
+
+  // Accumulate site-average plaquette (sum over planes / Nplanes / Nc)
+  LatticeComplex plaq_all(&grid);
+  plaq_all = Zero();
+  int Nplanes = 0;
+
+  // Per-site running sum over planes (Re Tr/Nc, normalised per plane)
+  LatticeRealD plaq_sum(&grid);
+  plaq_sum = Zero();
+
+  std::cout << GridLogMessage
+            << "======== Per-plane plaquette statistics ========" << std::endl;
+  std::cout << GridLogMessage
+            << std::setw(6)  << "plane"
+            << std::setw(20) << "max"
+            << std::setw(20) << "min"
+            << std::setw(20) << "average"
+            << std::endl;
+
+  for (int mu = 1; mu < Nd; mu++) {
+    for (int nu = 0; nu < mu; nu++) {
+
+      // Per-site trace of plaquette in (mu,nu) plane
+      LatticeComplex sitePlaq(&grid);
+      ColourWilsonLoops::traceDirPlaquette(sitePlaq, U, mu, nu);
+
+      plaq_all = plaq_all + sitePlaq;
+      Nplanes++;
+
+      // --- global average via sum() ---
+      TComplex Tsum  = sum(sitePlaq);
+      ComplexD csum  = TensorRemove(Tsum);
+      RealD    avg   = csum.real() / vol / Nc;
+
+      // --- global max and min via unvectorize + GlobalMax/GlobalMin ---
+      std::vector<TComplex> sv;
+      unvectorizeToLexOrdArray(sv, sitePlaq);
+
+      RealD local_max = -1e38, local_min = 1e38;
+      for (auto& tc : sv) {
+        RealD r = TensorRemove(tc).real() / Nc;
+        if (r > local_max) local_max = r;
+        if (r < local_min) local_min = r;
+      }
+
+      // Reduce across MPI ranks (no GlobalMin; negate to use GlobalMax)
+      grid.GlobalMax(local_max);
+      local_min = -local_min;
+      grid.GlobalMax(local_min);
+      local_min = -local_min;
+
+      std::cout << GridLogMessage
+                << std::setw(6)  << planeName(mu, nu)
+                << std::setw(20) << std::setprecision(10) << local_max
+                << std::setw(20) << std::setprecision(10) << local_min
+                << std::setw(20) << std::setprecision(10) << avg
+                << std::endl;
+
+      if (doThresh) {
+        std::vector<RealD> vr(sv.size(), 0.0);
+        for (int s = 0; s < (int)sv.size(); s++) {
+          RealD val = TensorRemove(sv[s]).real() / Nc;
+          if (val < threshold) {
+            Coordinate lc(Nd), gc(Nd);
+            Lexicographic::CoorFromIndex(lc, s, grid._ldimensions);
+            for (int d = 0; d < Nd; d++)
+              gc[d] = grid._processor_coor[d] * grid._ldimensions[d] + lc[d];
+            std::cerr << "BELOW_THRESHOLD  plane=" << planeName(mu, nu)
+//                      << "  global=(" << grid._processor_coor[0] << "," << grid._processor_coor[1] << "," << grid._processor_coor[2] << "," << grid._processor_coor[3] << ")"
+//                      << "  dim=(" << grid._ldimensions[0] << "," << grid._ldimensions[1] << "," << grid._ldimensions[2] << "," << grid._ldimensions[3] << ")"
+                      << "  site=(" << gc[0] << "," << gc[1] << "," << gc[2] << "," << gc[3] << ")"
+                      << "  P=" << std::setprecision(10) << val << std::endl;
+            vr[s] = val;
+          }
+        }
+        LatticeRealD plaqRe(&grid);
+        vectorizeFromLexOrdArray(vr, plaqRe);
+        plaq_sum = plaq_sum + plaqRe;
+      }
+    }
+  }
+
+  // Write per-site sum of below-threshold plaquette values over all planes
+  if (doThresh)
+    writeFile(plaq_sum, "./plaq_sum");
+
+  // ---- Overall (averaged over all planes) statistics ----
+  // plaq_all = sum of Tr[...] over all 6 (mu,nu) planes
+  // Normalise to Re Tr / Nc per plane
+  plaq_all = plaq_all * (1.0 / Nc / Nplanes);
+
+  TComplex Tsum_all = sum(plaq_all);
+  RealD    avg_all  = TensorRemove(Tsum_all).real() / vol;
+
+  std::vector<TComplex> sv_all;
+  unvectorizeToLexOrdArray(sv_all, plaq_all);
+
+  RealD max_all = -1e38, min_all = 1e38;
+  for (auto& tc : sv_all) {
+    RealD r = TensorRemove(tc).real();
+    if (r > max_all) max_all = r;
+    if (r < min_all) min_all = r;
+  }
+  grid.GlobalMax(max_all);
+  min_all = -min_all;
+  grid.GlobalMax(min_all);
+  min_all = -min_all;
+
+  // Cross-check with built-in avgPlaquette
+  RealD avg_builtin = ColourWilsonLoops::avgPlaquette(Umu);
+
+  std::cout << GridLogMessage
+            << "======== Overall plaquette statistics (all planes) ========" << std::endl;
+  std::cout << GridLogMessage << "  max     = " << max_all  << std::endl;
+  std::cout << GridLogMessage << "  min     = " << min_all  << std::endl;
+  std::cout << GridLogMessage << "  average = " << avg_all  << std::endl;
+  std::cout << GridLogMessage << "  avgPlaquette (builtin check) = " << avg_builtin << std::endl;
+
+  Grid_finalize();
+  return 0;
+}

tests/lanczos/Test_blocked_krylov_schur.cc

@@ -0,0 +1,99 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Test for BlockedKrylovSchur: verifies the KS decomposition A V = V H + F B^dag
+    by explicit operator applies, before and after each restart.
+
+    Uses DumbOperator (diagonal, real, Hermitian) from Test_synthetic_lanczos.
+
+    Tests Nblock=1 (scalar, regression) and Nblock=2,3 (exercises the B^H fix).
+
+*************************************************************************************/
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+// Diagonal real Hermitian operator (eigenvalues = scale lattice sites)
+template<class Field>
+class DumbOperator : public LinearOperatorBase<Field> {
+public:
+  LatticeComplex scale;
+
+  DumbOperator(GridBase* grid) : scale(grid) {
+    GridParallelRNG pRNG(grid);
+    std::vector<int> seeds({5,6,7,8});
+    pRNG.SeedFixedIntegers(seeds);
+    random(pRNG, scale);
+    scale = exp(-Grid::real(scale) * 3.0);
+  }
+
+  void OpDirAll(const Field& in, std::vector<Field>& out) {}
+  void OpDiag(const Field& in, Field& out) {}
+  void OpDir(const Field& in, Field& out, int dir, int disp) {}
+
+  void Op(const Field& in, Field& out) { out = scale * in; }
+  void AdjOp(const Field& in, Field& out) { out = scale * in; }
+  void HermOp(const Field& in, Field& out) { out = scale * in; }
+  void HermOpAndNorm(const Field& in, Field& out, double& n1, double& n2) {
+    out = scale * in;
+    ComplexD d = innerProduct(in, out);  n1 = real(d);
+    d = innerProduct(out, out);          n2 = real(d);
+  }
+};
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
+    GridDefaultLatt(),
+    GridDefaultSimd(Nd, vComplex::Nsimd()),
+    GridDefaultMpi());
+
+  GridParallelRNG RNG(grid);
+  RNG.SeedFixedIntegers({1,2,3,4});
+
+  typedef LatticeComplex Field;
+
+  DumbOperator<Field> op(grid);
+
+  int nFail = 0;
+
+  //--------------------------------------------------------------------
+  // Helper lambda: run BKS with doVerify and check it doesn't crash
+  //--------------------------------------------------------------------
+  auto runTest = [&](const std::string& label, int Nblock, int Nm, int Nk,
+                     int maxIter, int Nstop) {
+    std::cout << GridLogMessage << "===== " << label << " =====" << std::endl;
+
+    BlockKrylovSchur<Field> bks(op, grid, 1e-6, EvalReSmall);
+
+    std::vector<Field> v0(Nblock, Field(grid));
+    for (int t = 0; t < Nblock; t++) random(RNG, v0[t]);
+
+    bks(v0, maxIter, Nm, Nk, Nstop, Nblock,
+        /*doubleOrthog=*/true, /*doVerify=*/true);
+
+    std::cout << GridLogMessage << label << " done." << std::endl;
+  };
+
+  // Test 1: Nblock=1 — scalar case, regression (Nm,Nk now total vectors)
+  runTest("Nblock=1  Nm=10 Nk=5  maxIter=3", 1, 10, 5, 3, 5);
+
+  // Test 2: Nblock=2 — exercises the B^H fix for off-diagonal elements
+  runTest("Nblock=2  Nm=16 Nk=8  maxIter=3", 2, 16, 8, 3, 4);
+
+  // Test 3: Nblock=3 — further stress-test the B^H fix
+  runTest("Nblock=3  Nm=27 Nk=9  maxIter=3", 3, 27, 9, 3, 3);
+
+  // Test 4: Nblock=2, larger cycle — more restarts
+  runTest("Nblock=2  Nm=24 Nk=12  maxIter=5", 2, 24, 12, 5, 6);
+
+  if (nFail == 0)
+    std::cout << GridLogMessage << "All BlockedKrylovSchur tests completed." << std::endl;
+
+  Grid_finalize();
+  return nFail;
+}

tests/lanczos/Test_gamma5_block_lanczos.cc

@@ -0,0 +1,89 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Test for Gamma5BlockLanczos on a simple diagonal operator.
+
+    The operator is D = diag(scale_i) where scale is complex random.
+    γ5 is taken to be the identity (scalar field has no spin structure),
+    so the γ5-inner product reduces to the standard Euclidean inner product
+    and the algorithm should find the eigenvalues of D directly.
+
+    For a genuine Wilson Dirac test, pass the actual γ5 functor.
+
+*************************************************************************************/
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+// Diagonal complex operator: out = scale * in
+template<class Field>
+class DumbOperator : public LinearOperatorBase<Field> {
+public:
+  LatticeComplex scale;
+  DumbOperator(GridBase* grid) : scale(grid) {
+    GridParallelRNG pRNG(grid);
+    pRNG.SeedFixedIntegers({5,6,7,8});
+    random(pRNG, scale);
+    scale = exp(-Grid::real(scale) * 3.0);
+  }
+  void OpDirAll(const Field& in, std::vector<Field>& out) {}
+  void OpDiag(const Field& in, Field& out) {}
+  void OpDir(const Field& in, Field& out, int dir, int disp) {}
+  void Op(const Field& in, Field& out) { out = scale * in; }
+  void AdjOp(const Field& in, Field& out) { out = scale * in; }
+  void HermOp(const Field& in, Field& out) { out = scale * in; }
+  void HermOpAndNorm(const Field& in, Field& out, double& n1, double& n2) {
+    out = scale * in;
+    ComplexD d = innerProduct(in, out); n1 = real(d);
+    d = innerProduct(out, out);         n2 = real(d);
+  }
+};
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
+    GridDefaultLatt(),
+    GridDefaultSimd(Nd, vComplex::Nsimd()),
+    GridDefaultMpi());
+
+  GridParallelRNG RNG(grid);
+  RNG.SeedFixedIntegers({1,2,3,4});
+
+  typedef LatticeComplex Field;
+  DumbOperator<Field> op(grid);
+
+  // For LatticeComplex (scalar field) γ5 = identity
+  auto gamma5 = [](const Field& in, Field& out){ out = in; };
+
+  Field v0(grid);
+  random(RNG, v0);
+
+  const int maxSteps = 20;
+  const int Nstop    = 4;
+  const RealD tol    = 1e-6;
+
+  std::cout << GridLogMessage
+            << "\n========================================" << std::endl;
+  std::cout << GridLogMessage
+            << " Gamma5BlockLanczos  (maxSteps=" << maxSteps
+            << " Nstop=" << Nstop << ")" << std::endl;
+  std::cout << GridLogMessage
+            << "========================================\n" << std::endl;
+
+  Gamma5BlockLanczos<Field> g5bl(op, grid, gamma5, tol);
+  g5bl.doEvalCheck = true;
+  g5bl(v0, maxSteps, Nstop, /*reorthog=*/true);
+
+  auto evals = g5bl.getEvals();
+  std::cout << GridLogMessage
+            << "Gamma5BlockLanczos eigenvalues (" << evals.size() << "):" << std::endl;
+  for (int k = 0; k < (int)evals.size(); k++)
+    std::cout << GridLogMessage << "  [" << k << "] " << evals(k) << std::endl;
+
+  Grid_finalize();
+  return 0;
+}

tests/lanczos/Test_harmonic_vs_krylov_schur.cc

@@ -0,0 +1,154 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Comparison test: HarmonicBlockedKrylovSchur vs BlockedKrylovSchur.
+
+    Both algorithms are run on the same diagonal Hermitian operator and the
+    resulting eigenvalues are compared.  doVerify=true is used so the KS
+    decomposition check  max|H-M|  and the per-column residuals are printed
+    at each step.  For BKS these should be O(machine epsilon) at all times.
+    For HBKS they should be O(machine epsilon) AFTER Arnoldi, but may show
+    large per-column deviations AFTER restart+truncation (because the rotation
+    Q from Schur(Hhat) does not give an upper-triangular H_new, so the
+    truncated KS relation is only approximate).
+
+*************************************************************************************/
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+// Diagonal real Hermitian operator  out = scale * in
+template<class Field>
+class DumbOperator : public LinearOperatorBase<Field> {
+public:
+  LatticeComplex scale;
+  DumbOperator(GridBase* grid) : scale(grid) {
+    GridParallelRNG pRNG(grid);
+    pRNG.SeedFixedIntegers({5,6,7,8});
+    random(pRNG, scale);
+    scale = exp(-Grid::real(scale) * 3.0);
+  }
+  void OpDirAll(const Field& in, std::vector<Field>& out) {}
+  void OpDiag(const Field& in, Field& out) {}
+  void OpDir(const Field& in, Field& out, int dir, int disp) {}
+  void Op(const Field& in, Field& out) { out = scale * in; }
+  void AdjOp(const Field& in, Field& out) { out = scale * in; }
+  void HermOp(const Field& in, Field& out) { out = scale * in; }
+  void HermOpAndNorm(const Field& in, Field& out, double& n1, double& n2) {
+    out = scale * in;
+    ComplexD d = innerProduct(in, out); n1 = real(d);
+    d = innerProduct(out, out);         n2 = real(d);
+  }
+};
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
+    GridDefaultLatt(),
+    GridDefaultSimd(Nd, vComplex::Nsimd()),
+    GridDefaultMpi());
+
+  GridParallelRNG RNG(grid);
+  RNG.SeedFixedIntegers({1,2,3,4});
+
+  typedef LatticeComplex Field;
+  DumbOperator<Field> op(grid);
+
+  //----------------------------------------------------------------------
+  // Parameters (kept small so output is readable)
+  //----------------------------------------------------------------------
+  const int Nblock  = 4;
+  const int Nm      = 20;   // total vectors (= 5 blocks * Nblock=4)
+  const int Nk      = 8;    // total kept   (= 2 blocks * Nblock=4)
+  const int Nstop   = 4;
+  const int maxIter = 8;
+  const RealD tol   = 1e-6;
+
+  // Two identical starting blocks
+  std::vector<Field> v0(Nblock, Field(grid));
+  std::vector<Field> v0b(Nblock, Field(grid));
+  for (int t = 0; t < Nblock; t++) {
+    random(RNG, v0[t]);
+    v0b[t] = v0[t];   // identical start for fair comparison
+  }
+
+  //----------------------------------------------------------------------
+  // Run BlockedKrylovSchur with doVerify=true
+  //----------------------------------------------------------------------
+  std::cout << GridLogMessage
+            << "\n========================================" << std::endl;
+  std::cout << GridLogMessage
+            << " BlockKrylovSchur  (Nblock=" << Nblock
+            << " Nm=" << Nm << " Nk=" << Nk << ")" << std::endl;
+  std::cout << GridLogMessage
+            << "========================================\n" << std::endl;
+
+  BlockKrylovSchur<Field> bks(op, grid, tol, EvalImNormSmall);
+  bks(v0, maxIter, Nm, Nk, Nstop, Nblock,
+      /*doubleOrthog=*/true, /*doVerify=*/true);
+
+  auto bks_evals = bks.getEvals();
+  std::cout << GridLogMessage
+            << "BKS eigenvalues (" << bks_evals.size() << "):" << std::endl;
+  for (int k = 0; k < (int)bks_evals.size(); k++)
+    std::cout << GridLogMessage << "  [" << k << "] " << bks_evals[k] << std::endl;
+
+  //----------------------------------------------------------------------
+  // Run HarmonicBlockedKrylovSchur with doVerify=true
+  //----------------------------------------------------------------------
+  std::cout << GridLogMessage
+            << "\n========================================" << std::endl;
+  std::cout << GridLogMessage
+            << " HarmonicBlockKrylovSchur  (Nblock=" << Nblock
+            << " Nm=" << Nm << " Nk=" << Nk << " shift=0)" << std::endl;
+  std::cout << GridLogMessage
+            << "========================================\n" << std::endl;
+
+  HarmonicBlockKrylovSchur<Field> hbks(op, grid, tol, 0.0, EvalImNormSmall);
+  hbks(v0b, maxIter, Nm, Nk, Nstop, Nblock,
+       /*doubleOrthog=*/true, /*doVerify=*/true);
+
+  auto hbks_evals = hbks.getEvals();
+  std::cout << GridLogMessage
+            << "HBKS eigenvalues (" << hbks_evals.size() << "):" << std::endl;
+  for (int k = 0; k < (int)hbks_evals.size(); k++)
+    std::cout << GridLogMessage << "  [" << k << "] " << hbks_evals[k] << std::endl;
+
+  //----------------------------------------------------------------------
+  // Compare
+  //----------------------------------------------------------------------
+  std::cout << GridLogMessage
+            << "\n========================================" << std::endl;
+  std::cout << GridLogMessage << " Eigenvalue comparison" << std::endl;
+  std::cout << GridLogMessage
+            << "========================================" << std::endl;
+
+  // Sort both sets by real part for comparison
+  std::vector<ComplexD> bvec(bks_evals.data(),
+                             bks_evals.data() + bks_evals.size());
+  std::vector<ComplexD> hvec(hbks_evals.data(),
+                             hbks_evals.data() + hbks_evals.size());
+  auto cmpRe = [](const ComplexD& a, const ComplexD& b){ return a.real() < b.real(); };
+  std::sort(bvec.begin(), bvec.end(), cmpRe);
+  std::sort(hvec.begin(), hvec.end(), cmpRe);
+
+  int nCmp = std::min(bvec.size(), hvec.size());
+  double maxDiff = 0.0;
+  for (int k = 0; k < nCmp; k++) {
+    double diff = std::abs(bvec[k].real() - hvec[k].real()) + std::abs(bvec[k].imag() - hvec[k].imag());
+    maxDiff = std::max(maxDiff, diff);
+    std::cout << GridLogMessage
+              << "  k=" << k
+              << "  BKS=" << bvec[k]
+              << "  HBKS=" << hvec[k]
+              << "  |diff|=" << diff << std::endl;
+  }
+  std::cout << GridLogMessage << "  max |BKS - HBKS| = " << maxDiff << std::endl;
+
+  Grid_finalize();
+  return 0;
+}

tests/lanczos/Test_wilson_specflow.cc

@@ -200,8 +200,8 @@ int main(int argc, char** argv) {
     RD.close();
   }
 
-//  std::vector<Complex> boundary = {1,1,1,-1};
-  std::vector<Complex> boundary = {1,1,1,1};
+  std::vector<Complex> boundary = {1,1,1,-1};
+//  std::vector<Complex> boundary = {1,1,1,1};
   FermionOp::ImplParams Params(boundary);