Merge branch 'KS_shifted' of github.com:chulwoo1/Grid into KS_shifted

Grid/algorithms/iterative/BlockKrylovSchur.h

@@ -123,6 +123,7 @@ protected:
 
 public:
   std::vector<Field> evecs;
+  bool doEvalCheck = false;
 
   //--------------------------------------------------------------------
   // Constructor
@@ -252,6 +253,21 @@ public:
         std::cout << GridLogMessage << "BlockKrylovSchur: 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 << "BlockKrylovSchur: evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+
         return;
       }
     }
@@ -329,7 +345,7 @@ public:
     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::setw(14) << H(i, j);
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << H(i, j);
       std::cout << std::endl;
     }
 
@@ -337,7 +353,7 @@ public:
     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::setw(14) << B(i, t);
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << B(i, t);
       std::cout << std::endl;
     }
 
@@ -353,7 +369,7 @@ public:
     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::setw(14) << M(i, j);
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << M(i, j);
       std::cout << std::endl;
     }
 

Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h

@@ -34,65 +34,44 @@ See the full license in the file "LICENSE" in the top level distribution directo
 NAMESPACE_BEGIN(Grid);
 
 /**
- * Block harmonic restarted Krylov-Schur eigensolver.
+ * Block shift-targeted Krylov-Schur eigensolver.
  *
- * Harmonic Ritz values
+ * Algorithm
- * --------------------
+ * ---------
- * Standard Ritz values of A in a Krylov space K_m minimise the residual
+ * Uses a block Arnoldi factorisation:
- * in a Galerkin sense; they are good approximations to eigenvalues at the
- * *exterior* of the spectrum.  For eigenvalues *near* a target shift σ
- * (e.g. the smallest eigenvalues when σ=0) harmonic Ritz values are
- * better-suited: they are obtained by a Petrov-Galerkin condition that
- * requires the residual to be orthogonal to (A-σI)K_m instead of K_m.
- *
- * Given the block Arnoldi factorisation
  *
  *   A V = V H + F B^dag                                         (1)
  *
- * with V orthonormal (Nm columns), H the Nm² block
+ * with V orthonormal (Nm columns), H the Nm×Nm block upper-Hessenberg
- * upper-Hessenberg Rayleigh quotient, F the Nblock residual vectors and B
+ * Rayleigh quotient, F the Nblock residual vectors and B the Nm×Nblock
- * the Nm×Nblock coupling matrix, the harmonic Rayleigh quotient
+ * coupling matrix.
- * relative to shift σ is
- *
- *   Hhat = H + (H - σI)^{-H} B B^H                             (2)
- *
- * Derivation: the harmonic Ritz condition (A-σI)Vy ⊥ (A-σI)V leads to
- *
- *   [ (H-σI)^H (H-σI) + B B^H ] y = μ (H-σI)^H y
- *
- * Left-multiplying by (H-σI)^{-H} and setting θ = μ + σ gives the
- * standard eigenvalue problem  Hhat y = θ y  with Hhat as in (2).
- *
- * The harmonic Ritz values θ_j are eigenvalues of Hhat; among these,
- * the ones closest to σ (smallest |θ_j - σ|) are the best approximations
- * to the eigenvalues of A near σ.
  *
  * Thick restart
  * -------------
- * The Schur decomposition  Hhat = Q^dag S Q  is computed and the
+ * To target eigenvalues near shift σ, the Schur decomposition is computed
- * leading Nk*Nblock Schur values (sorted by the RitzFilter) are kept.
+ * for the shifted Rayleigh quotient:
- * The same unitary rotation Q is applied to both the Krylov basis and
- * to the *original* Rayleigh quotient H (not Hhat) for the restart:
  *
- *   V_new = V Q^dag [first Nk*Nblock columns]
+ *   (H - σI) = Q^dag S Q                                        (2)
- *   H_new = Q H Q^dag  [truncated Nk*Nblock × Nk*Nblock]
- *   B_new = Q B         [truncated Nk*Nblock × Nblock]
  *
- * Block Arnoldi then resumes from block Nk, restoring H to full size
+ * Sorting the Schur values of (H - σI) by smallest |S(i,i)| = |λ - σ| and
- * as new columns are appended.
+ * 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
  * -----------
- * For a harmonic Ritz pair (θ, y) the true Ritz residual bound is
+ * Convergence is declared when || B^H y_k || < Tolerance * approxLambdaMax
- *
+ * for each Ritz pair (λ_k, y_k) of the truncated H.
- *   || (A - θI) V y || ≤ || B^H y ||
- *
- * (same as for standard Ritz, because B captures the full coupling).
- * Convergence is declared when || B^H y || < Tolerance * approxLambdaMax.
  *
  * Parameters
  * ----------
- * shift    : target shift σ (default 0.0)
+ * 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)
@@ -146,6 +125,7 @@ class HarmonicBlockKrylovSchur {
 
 public:
   std::vector<Field>  evecs;
+  bool doEvalCheck = false;
 
   //--------------------------------------------------------------------
   // Constructor
@@ -205,33 +185,36 @@ public:
         verify(lbl);
       }
 
-      // ---- Form harmonic Rayleigh quotient ----
+      // ---- Schur decompose (H - σI) to select Schur vectors closest to σ ----
-      // Hhat = H + (H - σI)^{-H} * B * B^H
+      // Sorting the Schur values of (H - σI) by |S(i,i)| = |λ - σ| gives the
-      CMat Hhat = harmonicRayleigh(H, B, N);
+      // Ritz values of H nearest the target shift without any matrix inversion.
-
+      // Because Q diagonalises (H - σI), it also diagonalises H:
-      // ---- Schur decompose Hhat ----
+      //   H_new = Q H Q^dag = S + σI  (upper triangular)
-      ComplexSchurDecomposition schur(Hhat, false, ritzFilter);
+      // Truncation is therefore exact.
+      CMat Hshift = H - shift * CMat::Identity(N, N);
+      ComplexSchurDecomposition schur(Hshift, false, ritzFilter);
       schur.schurReorder(Nk);
 
       std::cout << GridLogMessage
-                << "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk):" << std::endl;
+                << "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) << std::endl;
+        std::cout << GridLogMessage << "  [" << i << "] " << S(i, i) + shift << std::endl;
 
       CMat Q  = schur.getMatrixQ();
       CMat Qt = Q.adjoint();
 
-      // ---- Rotate Krylov basis using Q from Hhat ----
+      // ---- Rotate Krylov basis ----
       std::vector<Field> basis2;
       constructUR(basis2, basis, Qt, N);
       basis = basis2;
 
-      // ---- Update H and B (rotate H, not Hhat) ----
+      // ---- Update H and B ----
-      H = Q * H * Qt;
+      // 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 ----
+      // ---- Truncate to Nk (exact: H upper triangular) ----
       int Nkeep = Nk;
 
       CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
@@ -249,9 +232,7 @@ public:
       std::cout << GridLogMessage
                 << "HarmonicBlockKrylovSchur: beta_k = " << beta_k << std::endl;
 
-      // Restart from the residual block F (unchanged from last Arnoldi step).
+      // Restart from F (exact: no discarded-basis correction needed)
-      // Note: for a Hermitian operator the correct H rows H[i,j] for i >= Nkeep+Nblock,
-      // j < Nkeep are filled via Hermitian symmetry inside blockArnoldiStep.
       startBlock = F;
 
       if (doVerify) {
@@ -273,6 +254,21 @@ public:
                   << "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;
       }
     }
@@ -292,8 +288,7 @@ public:
    *   A V = V H + F B^dag                                      (KS)
    *
    * by explicit operator applications.  H here is the standard Rayleigh
-   * quotient (not Hhat), so the KS relation is the same as for
+   * quotient, so the KS relation is the same as for BlockKrylovSchur.
-   * BlockKrylovSchur.
    *
    * Prints:
    *   - H (current Rayleigh quotient, nBasis × nBasis)
@@ -326,7 +321,7 @@ public:
     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::setw(14) << H(i, j);
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << H(i, j);
       std::cout << std::endl;
     }
 
@@ -334,7 +329,7 @@ public:
     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::setw(14) << B(i, t);
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << B(i, t);
       std::cout << std::endl;
     }
 
@@ -350,7 +345,7 @@ public:
     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::setw(14) << M(i, j);
+        std::cout << "  " << std::setprecision(4) << std::setw(14) << M(i, j);
       std::cout << std::endl;
     }
 
@@ -410,32 +405,6 @@ public:
 
 private:
 
-  //--------------------------------------------------------------------
-  // Harmonic Rayleigh quotient
-  //--------------------------------------------------------------------
-  /**
-   * Forms the harmonic Rayleigh quotient relative to shift σ:
-   *
-   *   Hhat = H + (H - σI)^{-H} * B * B^H
-   *
-   * where H is the N×N block-Hessenberg, B is the N×Nblock coupling matrix.
-   *
-   * The N×N solve  (H - σI)^H X = B B^H  is done via Eigen's LU
-   * factorisation.  If H - σI is (nearly) singular the result is
-   * ill-conditioned; in that case σ should be perturbed slightly.
-   */
-  CMat harmonicRayleigh(const CMat& H_, const CMat& B_, int N)
-  {
-    CMat K  = H_ - shift * CMat::Identity(N, N);
-    CMat KH = K.adjoint();                         // (H - σI)^H
-
-    // Solve KH * X = B B^H  →  X = KH^{-1} B B^H = (H-σI)^{-H} B B^H
-    CMat BBH = B_ * B_.adjoint();                  // N × N
-    CMat X   = KH.lu().solve(BBH);                 // N × N
-
-    return H_ + X;
-  }
-
   //--------------------------------------------------------------------
   // Block Arnoldi iteration
   //--------------------------------------------------------------------
@@ -454,25 +423,17 @@ private:
       blockOrthonormalise(V0);
       for (auto& v : V0) basis.push_back(v);
     } else {
-      // Append residual block (startBlock = F_old) to basis.
+      // Append the new starting block to the retained basis.
-      // The truncated KS relation after restart is:
       //
-      //   A V_k = V_k S_k + F_old B_old^dag          (*)
+      // 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]).
-      // where V_k = basis[0:Nkeep], S_k is stored in H[0:Nkeep,0:Nkeep],
+      // Since H_new = S + σI is upper triangular, the off-diagonal block
-      // B_old = B[0:Nkeep,:], F_old = startBlock.
+      // H_new[Nkeep:, :Nkeep] = 0 and the restart from F is exact.
-      //
-      // 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 "restart coupling rows" must be set before Arnoldi continues.
       int Nkeep = startIdx * Nblock;
       for (auto& v : startBlock) basis.push_back(v);
 
-      // Fill restart coupling rows into H
+      // 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));
@@ -524,10 +485,6 @@ private:
         for (int s = 0; s < Nblock; s++)
           B(kBase + t, s) = std::conj(R(s, t));    // B_block = R^H
       beta_k = R.norm();
-      // Hermitian symmetry fill for last block (same as non-last path below)
-      for (int t = 0; t < Nblock; t++)
-        for (int j = 0; j < kBase; j++)
-          H(kBase + t, j) = std::conj(H(j, kBase + t));
       return;
     }
 
@@ -542,23 +499,6 @@ private:
     for (int t = 0; t < Nblock; t++)
       basis.push_back(F[t]);
 
-    // Hermitian symmetry fill: H[kBase+t, j] = conj(H[j, kBase+t]) for j < kBase.
-    //
-    // In a fresh block Arnoldi the Krylov structure forces H[kBase+t, j] = 0 for
-    // j < kBase-Nblock (sub-subdiagonal), so this is a no-op.
-    //
-    // After a non-Schur restart (e.g. harmonic restart where H_new = Q H Q^dag is
-    // a full matrix), A v_k_j for j < Nkeep has components in ALL new extended
-    // vectors, making these elements non-zero.  The Arnoldi step fills column
-    // kBase+t (H[j, kBase+t] for j < prevN) via inner products, but never fills
-    // the corresponding row.  For a Hermitian operator the two are related by
-    //   H[kBase+t, j] = <basis[kBase+t] | A basis[j]>
-    //                 = conj(<basis[j] | A basis[kBase+t]>) = conj(H[j, kBase+t])
-    // Filling these ensures H = H^dag and fixes the M != H discrepancy that
-    // corrupts subsequent Arnoldi steps after a harmonic restart.
-    for (int t = 0; t < Nblock; t++)
-      for (int j = 0; j < kBase; j++)
-        H(kBase + t, j) = std::conj(H(j, kBase + t));
   }
 
   //--------------------------------------------------------------------

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
@@ -338,6 +339,7 @@ class KrylovSchur {
 
   public:
     RealD *shift;
+    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),
@@ -604,7 +606,21 @@ if (!shift){
           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;
+          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;
+            }
+          }
 
           // writeEigensystem(path);
 

examples/Example_krylov_schur.cc

@@ -56,6 +56,7 @@ struct LanczosParameters: Serializable {
                                 Integer, maxIter,
                                 Integer, Nblock,
                                 Integer, verify,
+		  		RealD, shift , 
 	  			RealD, resid,
 	  			RealD, ChebyLow,
 	  			RealD, ChebyHigh,
@@ -117,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){
@@ -131,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);
@@ -340,29 +333,22 @@ int main (int argc, char ** argv)
   LatticeFermion src2 = src[0];
 
   // Run KrylovSchur and Arnoldi on a Hermitian matrix
-  std::cout << GridLogMessage << "Running Krylov Schur" << std::endl;
+    RealD shift=LanParams.shift;
-  // 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
-    RealD shift=1.5;
 #if 1
+    std::cout << GridLogMessage << "Running Krylov Schur" << std::endl;
     KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
 //    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
     KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
     std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
 #else
+    std::cout << GridLogMessage << "Running BlockKrylovSchur" << std::endl;
     int Nblock=4;
     Nblock=LanParams.Nblock;
     bool if_verify=false;
     if(LanParams.verify) if_verify=true;
-//    KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+//    BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
-//    KrySchur(src, maxIter, Nm, Nk, Nstop,true,if_verify);
+    HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall);
-    BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
     KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify);
-//    HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall);
-//    KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true);
     std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl;
 #endif
 

tests/core/Test_plaquette_stats.cc

@@ -0,0 +1,215 @@
+/*************************************************************************************
+
+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]
+ *
+ *   --file <path>  Read gauge field from NERSC-format file
+ *   --hot          Use random (hot) SU(3) start (default: cold/unit start)
+ *
+ * Grid size defaults to 8^3 x 16; override with --grid (e.g. --grid 4.4.4.8).
+ */
+
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+// 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; }
+  }
+
+  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;
+
+  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;
+    }
+  }
+
+  // ---- 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_harmonic_vs_krylov_schur.cc

@@ -61,11 +61,11 @@ int main(int argc, char** argv)
   //----------------------------------------------------------------------
   // Parameters (kept small so output is readable)
   //----------------------------------------------------------------------
-  const int Nblock  = 2;
+  const int Nblock  = 4;
-  const int Nm      = 12;   // total vectors (= 6 blocks * Nblock=2)
+  const int Nm      = 20;   // total vectors (= 5 blocks * Nblock=4)
-  const int Nk      = 6;    // total kept   (= 3 blocks * Nblock=2)
+  const int Nk      = 8;    // total kept   (= 2 blocks * Nblock=4)
-  const int Nstop   = 2;
+  const int Nstop   = 4;
-  const int maxIter = 4;
+  const int maxIter = 8;
   const RealD tol   = 1e-6;
 
   // Two identical starting blocks
@@ -73,7 +73,7 @@ int main(int argc, char** argv)
   std::vector<Field> v0b(Nblock, Field(grid));
   for (int t = 0; t < Nblock; t++) {
     random(RNG, v0[t]);
-    v0b[t] = v0[t];
+    v0b[t] = v0[t];   // identical start for fair comparison
   }
 
   //----------------------------------------------------------------------
@@ -87,7 +87,7 @@ int main(int argc, char** argv)
   std::cout << GridLogMessage
             << "========================================\n" << std::endl;
 
-  BlockKrylovSchur<Field> bks(op, grid, tol, EvalReSmall);
+  BlockKrylovSchur<Field> bks(op, grid, tol, EvalImNormSmall);
   bks(v0, maxIter, Nm, Nk, Nstop, Nblock,
       /*doubleOrthog=*/true, /*doVerify=*/true);
 
@@ -108,7 +108,7 @@ int main(int argc, char** argv)
   std::cout << GridLogMessage
             << "========================================\n" << std::endl;
 
-  HarmonicBlockKrylovSchur<Field> hbks(op, grid, tol, 0.0, EvalNormSmall);
+  HarmonicBlockKrylovSchur<Field> hbks(op, grid, tol, 0.0, EvalImNormSmall);
   hbks(v0b, maxIter, Nm, Nk, Nstop, Nblock,
        /*doubleOrthog=*/true, /*doVerify=*/true);