From e8e7ef08fcb94ff092f562f8e52c36367b3129b2 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 15 Jan 2026 03:20:01 +0000 Subject: [PATCH] KrylovSchur and spectral flow updates --- examples/Example_krylov_schur.cc | 12 ++-- tests/lanczos/Test_wilson_specflow.cc | 86 ++++++++++++++++++++++++--- 2 files changed, 87 insertions(+), 11 deletions(-) diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index d194d13d..9d5d4d22 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -126,12 +126,16 @@ public: assert(0); }; 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 HermOp(_Mat); - BiCGSTAB CG(_stp,10000); - CG(HermOp,in,out); +// NonHermitianLinearOperator HermOp(_Mat); +// BiCGSTAB CG(_stp,10000); + _Mat.Mdag(in,tmp); + MdagMLinearOperator HermOp(_Mat); + ConjugateGradient CG(_stp,10000); + CG(HermOp,tmp,out); // out = out + shift * in; } void AdjOp (const Field &in, Field &out){ @@ -331,7 +335,7 @@ int main (int argc, char ** argv) // KrySchur(src, maxIter, Nm, Nk, Nstop); // KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); // Hacked, really EvalImagSmall -#if 0 +#if 1 RealD shift=1.5; KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); diff --git a/tests/lanczos/Test_wilson_specflow.cc b/tests/lanczos/Test_wilson_specflow.cc index e9bd04df..4fd63a70 100644 --- a/tests/lanczos/Test_wilson_specflow.cc +++ b/tests/lanczos/Test_wilson_specflow.cc @@ -27,6 +27,7 @@ directory *************************************************************************************/ /* END LEGAL */ #include +#include using namespace std; using namespace Grid; @@ -38,11 +39,32 @@ typedef typename WilsonFermionD::FermionField FermionField; RealD AllZero(RealD x) { return 0.; } +template void writeFile(T& in, std::string const fname){ +#if 1 + // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter WR(in.Grid()->IsBoss()); + WR.open(fname); + WR.writeScidacFieldRecord(in,record,0); + WR.close(); +#endif + // What is the appropriate way to throw error? +} + + namespace Grid { struct LanczosParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters, RealD, mass , + RealD, mstep , + Integer, Nstop, + Integer, Nk, + Integer, Np, + Integer, ReadEvec, + Integer, maxIter, + RealD, resid, RealD, ChebyLow, RealD, ChebyHigh, Integer, ChebyOrder) @@ -115,12 +137,13 @@ int main(int argc, char** argv) { LatticeGaugeField Umu(UGrid); // SU::HotConfiguration(RNG4, Umu); +// SU::ColdConfiguration(Umu); FieldMetaData header; std::string file("./config"); - int precision32 = 0; - int tworow = 0; +// int precision32 = 0; +// int tworow = 0; // NerscIO::writeConfiguration(Umu,file,tworow,precision32); NerscIO::readConfiguration(Umu,header,file); @@ -158,10 +181,33 @@ int main(int argc, char** argv) { } mass=LanParams.mass; + resid=LanParams.resid; + Nstop=LanParams.Nstop; + Nk=LanParams.Nk; + Np=LanParams.Np; + MaxIt=LanParams.maxIter; + Nm = Nk + Np; + + FermionField src(FGrid); + gaussian(RNG5, src); + if(LanParams.ReadEvec) { + std::string evecs_file="evec_in"; + std::cout << GridLogIRL<< "Reading evecs from "< boundary = {1,1,1,-1}; +// std::vector boundary = {1,1,1,0}; +// std::vector boundary = {1,1,1,1}; + FermionOp::ImplParams Params(boundary); -while ( mass > - 5.0){ - FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass); +while ( mass > - 2.0){ + FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass,Params); MdagMLinearOperator HermOp(WilsonOperator); /// <----- //SchurDiagTwoOperator HermOp(WilsonOperator); Gamma5HermitianLinearOperator HermOp2(WilsonOperator); /// <----- @@ -180,10 +226,9 @@ while ( mass > - 5.0){ PlainHermOp Op2 (HermOp2); ImplicitlyRestartedLanczos IRL(OpCheby, Op2, Nstop, Nk, Nm, resid, MaxIt); +// SimpleLanczos IRL(Op,Nstop, Nk, Nm, resid, MaxIt); std::vector eval(Nm); - FermionField src(FGrid); - gaussian(RNG5, src); std::vector evec(Nm, FGrid); for (int i = 0; i < 1; i++) { std::cout << i << " / " << Nm << " grid pointer " << evec[i].Grid() @@ -192,19 +237,46 @@ while ( mass > - 5.0){ int Nconv; IRL.calc(eval, evec, src, Nconv); +// IRL.calc(eval, src, Nconv); std::cout << mass <<" : " << eval << std::endl; Gamma g5(Gamma::Algebra::Gamma5) ; ComplexD dot; FermionField tmp(FGrid); + FermionField sav(FGrid); + sav=evec[0]; for (int i = 0; i < Nstop ; i++) { tmp = g5*evec[i]; dot = innerProduct(tmp,evec[i]); std::cout << mass << " : " << eval[i] << " " << real(dot) << " " << imag(dot) << std::endl ; +// if ( i<1) + { + std::string evfile ("./evec_"+std::to_string(mass)+"_"+std::to_string(i)); + auto evdensity = localInnerProduct(evec[i],evec[i] ); + writeFile(evdensity,evfile); +// if(LanParams.ReadEvec) { +// std::string evecs_file="evec_in"; + { + std::cout << GridLogIRL<< "Reading evecs from "<0) sav += evec[i]; + } + { + std::string evfile ("./evec_"+std::to_string(mass)+"_sum"); +// auto evdensity = localInnerProduct(evec[i],evec[i] ); + writeFile(sav,evfile); } src = evec[0]+evec[1]+evec[2]; - mass += -0.1; + src += evec[3]+evec[4]+evec[5]; + src += evec[6]+evec[7]+evec[8]; + mass += LanParams.mstep; } Grid_finalize();