From 786496f22e91cf9299b67c1566d3344618a1e1bc Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Mon, 3 Nov 2025 21:18:56 +0000 Subject: [PATCH] Checking in before pulling KrylovSchur --- .../iterative/ImplicitlyRestartedLanczos.h | 27 +- tests/lanczos/Test_wilson_lanczos.cc | 242 ++++++++++++++++-- tests/lanczos/Test_wilson_specflow.cc | 37 ++- 3 files changed, 275 insertions(+), 31 deletions(-) diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h index eb0d0761..df2007d2 100644 --- a/Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h +++ b/Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h @@ -53,6 +53,18 @@ enum IRLdiagonalisation { IRLdiagonaliseWithEigen }; +enum IRLeigsort { + IRLeigsortMax, + IRLeigsortSqMin +}; + +#if 0 +bool square_comp(RealD a, RealD b){ + if (a*a class ImplicitlyRestartedLanczosHermOpTester : public ImplicitlyRestartedLanczosTester { public: @@ -119,9 +131,10 @@ class ImplicitlyRestartedLanczos { ///////////////////////// // Constructor ///////////////////////// - -public: + public: + IRLeigsort EigSort; + ////////////////////////////////////////////////////////////////// // PAB: ////////////////////////////////////////////////////////////////// @@ -154,6 +167,7 @@ public: Nstop(_Nstop) , Nk(_Nk), Nm(_Nm), eresid(_eresid), betastp(_betastp), MaxIter(_MaxIter) , MinRestart(_MinRestart), + EigSort(IRLeigsortMax), orth_period(_orth_period), diagonalisation(_diagonalisation) { }; ImplicitlyRestartedLanczos(LinearFunction & PolyOp, @@ -170,6 +184,7 @@ public: Nstop(_Nstop) , Nk(_Nk), Nm(_Nm), eresid(_eresid), betastp(_betastp), MaxIter(_MaxIter) , MinRestart(_MinRestart), + EigSort(IRLeigsortMax), orth_period(_orth_period), diagonalisation(_diagonalisation) { }; //////////////////////////////// @@ -316,8 +331,12 @@ until convergence // sorting ////////////////////////////////// eval2_copy = eval2; +// if (EigSort==IRLeigsortMax) +// std::partial_sort(eval2.begin(),eval2.begin()+Nm,eval2.end(),square_comp); +// else std::partial_sort(eval2.begin(),eval2.begin()+Nm,eval2.end(),std::greater()); std::cout< +Author: Chulwoo Jung 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 @@ -27,6 +27,9 @@ directory *************************************************************************************/ /* END LEGAL */ #include +#include +#include + using namespace std; using namespace Grid; @@ -38,18 +41,111 @@ 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, + RealD, resid, + RealD, ChebyLow, + RealD, ChebyHigh, + Integer, ChebyOrder) +// Integer, StartTrajectory, +// Integer, Trajectories, /* @brief Number of sweeps in this run */ +// bool, MetropolisTest, +// Integer, NoMetropolisUntil, +// std::string, StartingType, +// Integer, SW, +// RealD, Kappa, +// IntegratorParameters, MD) + + LanczosParameters() { + ////////////////////////////// Default values + mass = 0; +// MetropolisTest = true; +// NoMetropolisUntil = 10; +// StartTrajectory = 0; +// SW = 2; +// Trajectories = 10; +// StartingType = "HotStart"; + ///////////////////////////////// + } + + template + LanczosParameters(Reader & TheReader){ + initialize(TheReader); + } + + template < class ReaderClass > + void initialize(Reader &TheReader){ +// std::cout << GridLogMessage << "Reading HMC\n"; + read(TheReader, "HMC", *this); + } + + + void print_parameters() const { +// std::cout << GridLogMessage << "[HMC parameters] Trajectories : " << Trajectories << "\n"; +// std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n"; +// std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n"; +// std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n"; +// std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n"; +// MD.print_parameters(); + } + +}; + +} + int main(int argc, char** argv) { Grid_init(&argc, &argv); + int Ndir=4; + auto mpi_layout = GridDefaultMpi(); + std::vector nblock(4,1); + std::vector mpi_split(4,1); +//Interested in avoiding degeneracy only for now + nblock[3]=2; + + int mrhs=1; + for(int i =0;i seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); @@ -62,12 +158,15 @@ int main(int argc, char** argv) { LatticeGaugeField Umu(UGrid); // SU::HotConfiguration(RNG4, Umu); - SU::ColdConfiguration(Umu); +// SU::ColdConfiguration(Umu); FieldMetaData header; std::string file("./config"); - NerscIO::readConfiguration(Umu,header,file); +// int precision32 = 0; +// int tworow = 0; +// NerscIO::writeConfiguration(Umu,file,tworow,precision32); + NerscIO::readConfiguration(Umu,header,file); /* std::vector U(4, UGrid); @@ -75,38 +174,101 @@ int main(int argc, char** argv) { U[mu] = PeekIndex(Umu, mu); } */ + + int Nstop = 10; + int Nu = 1; + int Nk = 20; + int Np = 80; + int Nm = Nk + Np; + int MaxIt = 10000; + RealD resid = 1.0e-5; + + RealD mass = -1.0; + + LanczosParameters LanParams; +#if 1 + { + XmlReader HMCrd("LanParams.xml"); + read(HMCrd,"LanczosParameters",LanParams); + } +#else + { + LanParams.mass = mass; + } +#endif + std::cout << GridLogMessage<< LanParams < src(Nu,FGrid); + for(int i =0;i boundary = {1,1,1,-1}; // std::vector boundary = {1,1,1,1}; FermionOp::ImplParams Params(boundary); + GridCartesian * SFGrid = SGrid; + GridRedBlackCartesian * SFrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(SFGrid); +// GridRedBlackCartesian * SFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(JP.Ls,SGrid); + + LatticeGaugeField s_Umu(SGrid); + Grid_split (Umu,s_Umu); - RealD mass = 0.0; -// FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass); + +while ( mass > - 2.0){ FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass,Params); - MdagMLinearOperator HermOp(WilsonOperator); /// <----- + MdagMLinearOperator HermOp(WilsonOperator); /// <----- + FermionOp WilsonSplit(s_Umu,*SFGrid,*SFrbGrid,mass,Params); + MdagMLinearOperator SHermOp(WilsonSplit); /// <----- //SchurDiagTwoOperator HermOp(WilsonOperator); - - const int Nstop = 20; - const int Nk = 60; - const int Np = 60; - const int Nm = Nk + Np; - const int MaxIt = 10000; - RealD resid = 1.0e-6; + Gamma5HermitianLinearOperator HermOp2(WilsonOperator); /// <----- std::vector Coeffs{0, 1.}; Polynomial PolyX(Coeffs); - Chebyshev Cheby(0.0, 10., 12); +// Chebyshev Cheby(0.5, 60., 31); +// RealD, ChebyLow, +// RealD, ChebyHigh, +// Integer, ChebyOrder) + + Chebyshev Cheby(LanParams.ChebyLow,LanParams.ChebyHigh,LanParams.ChebyOrder); FunctionHermOp OpCheby(Cheby,HermOp); PlainHermOp Op (HermOp); + PlainHermOp Op2 (HermOp2); -// ImplicitlyRestartedLanczos IRL(OpCheby, Op, Nstop, Nk, Nm, resid, MaxIt); - SimpleLanczos IRL(Op,Nstop, Nk, Nm, resid, MaxIt); - +// ImplicitlyRestartedLanczos IRL(OpCheby, Op2, Nstop, Nk, Nm, resid, MaxIt); +// SimpleLanczos IRL(Op,Nstop, Nk, Nm, resid, MaxIt); + ImplicitlyRestartedBlockLanczos IRBL(HermOp, SHermOp, + FrbGrid,SFrbGrid,mrhs, + Cheby, + Nstop, Nstop*2, + Nu, Nk, Nm, + resid, MaxIt, + IRBLdiagonaliseWithEigen); + IRBL.split_test=1; 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() @@ -115,9 +277,39 @@ int main(int argc, char** argv) { int Nconv; // IRL.calc(eval, evec, src, Nconv); - IRL.calc(eval, src, Nconv); + IRBL.calc(eval, evec, src, Nconv,LanczosType::irbl); - std::cout << eval << std::endl; + 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 (i>0) sav += evec[i]; + } + { + std::string evfile ("./evec_"+std::to_string(mass)+"_sum"); +// auto evdensity = localInnerProduct(evec[i],evec[i] ); + writeFile(sav,evfile); + } + for(int i =0;i boundary = {1,1,1,-1}; // std::vector boundary = {1,1,1,1}; FermionOp::ImplParams Params(boundary); -while ( mass > - 2.5){ +while ( mass > - 2.0){ FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass,Params); MdagMLinearOperator HermOp(WilsonOperator); /// <----- //SchurDiagTwoOperator HermOp(WilsonOperator); @@ -228,6 +241,8 @@ while ( mass > - 2.5){ 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]); @@ -237,7 +252,23 @@ while ( mass > - 2.5){ 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]; src += evec[3]+evec[4]+evec[5];