From 3fac61fc558ae93397e889838e59c2900da1262b Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Fri, 29 May 2026 13:55:27 -0400
Subject: [PATCH] benchmarks: add DhopEO benchmarks to Benchmark_staggered and
 Benchmark_staggeredF

Also add NaiveStaggeredFermionF to Benchmark_staggeredF, and fix bug where
NaiveStaggered timed loop was calling Ds.Dhop instead of Dn.Dhop.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 benchmarks/Benchmark_staggered.cc  | 64 +++++++++++++++++-------------
 benchmarks/Benchmark_staggeredF.cc | 49 +++++++++++++++++------
 2 files changed, 72 insertions(+), 41 deletions(-)

diff --git a/benchmarks/Benchmark_staggered.cc b/benchmarks/Benchmark_staggered.cc
index a419d600..a9f6bc57 100644
--- a/benchmarks/Benchmark_staggered.cc
+++ b/benchmarks/Benchmark_staggered.cc
@@ -97,37 +97,45 @@ int main (int argc, char ** argv)
   RealD u0=1.0;
   ImprovedStaggeredFermionD Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
   NaiveStaggeredFermionD Dn(Umu,Grid,RBGrid,mass,c1,u0,params);
-  
-  std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
-  int ncall=100;
-  // warm perf only
-  for(int i=0;i<ncall;i++){
-    Ds.Dhop(src,result,0);
-  }
-  double t0=usecond();
-  for(int i=0;i<ncall;i++){
-    Ds.Dhop(src,result,0);
-  }
-  double t1=usecond();
-  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 + 90 == 1146
-  std::cout<<GridLogMessage << "Called Ds"<<std::endl;
-  std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
-  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
 
-  // Warm perf only
-  for(int i=0;i<ncall;i++){
-    Dn.Dhop(src,result,0);
-  }
+  FermionField src_e(&RBGrid); pickCheckerboard(Even, src_e, src);
+  FermionField src_o(&RBGrid); pickCheckerboard(Odd,  src_o, src);
+  FermionField res_e(&RBGrid); res_e = Zero();
+
+  int ncall=100;
+
+  // ImprovedStaggered Dhop
+  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
+  double t0=usecond();
+  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
+  double t1=usecond();
+  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // 1146 flops/site
+  std::cout<<GridLogMessage << "ImprovedStaggered Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
+  // ImprovedStaggered DhopEO
+  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
   t0=usecond();
-  for(int i=0;i<ncall;i++){
-    Ds.Dhop(src,result,0);
-  }
+  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
   t1=usecond();
-  flops=(8*(3*(6+8+8)) + 7*3*2)*volume*ncall; 
-  std::cout<<GridLogMessage << "Called Dn"<<std::endl;
-  std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
-  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
-  
+  flops=(16*(3*(6+8+8)) + 15*3*2)*(volume/2)*ncall;
+  std::cout<<GridLogMessage << "ImprovedStaggered DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
+  // NaiveStaggered Dhop
+  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
+  t0=usecond();
+  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
+  t1=usecond();
+  flops=(8*(3*(6+8+8)) + 7*3*2)*volume*ncall;
+  std::cout<<GridLogMessage << "NaiveStaggered    Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
+  // NaiveStaggered DhopEO
+  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
+  t0=usecond();
+  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
+  t1=usecond();
+  flops=(8*(3*(6+8+8)) + 7*3*2)*(volume/2)*ncall;
+  std::cout<<GridLogMessage << "NaiveStaggered    DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
   Grid_finalize();
 }
 
diff --git a/benchmarks/Benchmark_staggeredF.cc b/benchmarks/Benchmark_staggeredF.cc
index e0f4331f..a858be56 100644
--- a/benchmarks/Benchmark_staggeredF.cc
+++ b/benchmarks/Benchmark_staggeredF.cc
@@ -96,22 +96,45 @@ int main (int argc, char ** argv)
   RealD c2=-1.0/24.0;
   RealD u0=1.0;
   ImprovedStaggeredFermionF Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
-  
-  std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
+  NaiveStaggeredFermionF    Dn(Umu,Grid,RBGrid,mass,c1,u0,params);
+
+  FermionField src_e(&RBGrid); pickCheckerboard(Even, src_e, src);
+  FermionField src_o(&RBGrid); pickCheckerboard(Odd,  src_o, src);
+  FermionField res_e(&RBGrid); res_e = Zero();
+
   int ncall=1000;
-  for(int i=0;i<ncall;i++){
-    Ds.Dhop(src,result,0);
-  }
+
+  // ImprovedStaggered Dhop
+  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
   double t0=usecond();
-  for(int i=0;i<ncall;i++){
-    Ds.Dhop(src,result,0);
-  }
+  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
   double t1=usecond();
-  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 +  == 1146
-  
-  std::cout<<GridLogMessage << "Called Ds"<<std::endl;
-  std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
-  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
+  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // 1146 flops/site
+  std::cout<<GridLogMessage << "ImprovedStaggered Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
+  // ImprovedStaggered DhopEO
+  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
+  t0=usecond();
+  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
+  t1=usecond();
+  flops=(16*(3*(6+8+8)) + 15*3*2)*(volume/2)*ncall;
+  std::cout<<GridLogMessage << "ImprovedStaggered DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
+  // NaiveStaggered Dhop
+  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
+  t0=usecond();
+  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
+  t1=usecond();
+  flops=(8*(3*(6+8+8)) + 7*3*2)*volume*ncall;
+  std::cout<<GridLogMessage << "NaiveStaggered    Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
+
+  // NaiveStaggered DhopEO
+  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
+  t0=usecond();
+  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
+  t1=usecond();
+  flops=(8*(3*(6+8+8)) + 7*3*2)*(volume/2)*ncall;
+  std::cout<<GridLogMessage << "NaiveStaggered    DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
 
   Grid_finalize();
 }