From caa66418bd795826efda42812262c59d5da61379 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 6 Nov 2025 22:44:05 +0000 Subject: [PATCH 1/2] Checking in before pulling --- Grid/algorithms/iterative/Arnoldi.h | 13 +- Grid/algorithms/iterative/KrylovSchur.h | 25 ++- examples/Example_krylov_schur.cc | 197 ++++++++++++++++++++---- examples/Example_spec_kryschur.cc | 45 ++++++ 4 files changed, 240 insertions(+), 40 deletions(-) diff --git a/Grid/algorithms/iterative/Arnoldi.h b/Grid/algorithms/iterative/Arnoldi.h index 19b5f6c7..b4a9b8b8 100644 --- a/Grid/algorithms/iterative/Arnoldi.h +++ b/Grid/algorithms/iterative/Arnoldi.h @@ -51,14 +51,23 @@ struct ComplexComparator ComplexComparator (RitzFilter _f) : f(_f) {} bool operator()(std::complex z1, std::complex z2) { switch (f) { + RealD tmp1, tmp2; + tmp1=std::abs(std::imag(z1)); + tmp2=std::abs(std::imag(z2)); case EvalNormSmall: return std::abs(z1) < std::abs(z2); case EvalNormLarge: return std::abs(z1) > std::abs(z2); +// Terrible hack +// return std::abs(std::real(z1)) < std::abs(std::real(z2)); +// if ( std::abs(std::real(z1)) >4.) tmp1 +=1.; +// if ( std::abs(std::real(z2)) >4.) tmp2 +=1.; case EvalReSmall: - return std::abs(std::real(z1)) < std::abs(std::real(z2)); + return tmp1 < tmp2; +// return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); case EvalReLarge: - return std::abs(std::real(z1)) > std::abs(std::real(z2)); + return tmp1 > tmp2; +// return std::abs(std::real(z1)) > std::abs(std::real(z2)); default: assert(0); } diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 0c660585..5441af1a 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -145,8 +145,8 @@ class ComplexSchurDecomposition { * TODO: pass in compare function as an argument, default to compare with <. */ // void schurReorder(int Nk, std::function compare) { - void schurReorder(int Nk) { - for (int i = 0; i < Nk; i++) { + int schurReorder(int Nk) { + for (int i = 0; i < Nm; i++) { for (int k = 0; k <= Nm - 2; k++) { int idx = Nm - 2 - k; // TODO use RitzFilter enum here @@ -157,7 +157,12 @@ class ComplexSchurDecomposition { } } } - return; + int temp=1; + for (int i = 0; i < Nm-1; i++) { + std::cout << GridLogMessage << "S " << i << " "< evecs; // Vector of evec fields RitzFilter ritzFilter; // how to sort evals @@ -210,6 +214,8 @@ class KrylovSchur { u = Zero(); }; + std::vector evecs; // Vector of evec fields + /* Getters */ Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; } @@ -257,8 +263,10 @@ class KrylovSchur { std::cout << GridLogDebug << "Schur decomp holds? " << schur.checkDecomposition() << std::endl; // Rearrange Schur matrix so wanted evals are on top left (like MATLAB's ordschur) - std::cout << GridLogMessage << "Reordering Schur eigenvalues" << std::endl; + std::cout << GridLogMessage << "Reordering Schur eigenvalues Nk=" << Nk << std::endl; +// Nk=schur.schurReorder(_Nk); schur.schurReorder(Nk); + std::cout << GridLogMessage << "After Reorder Nk= "< using namespace std; using namespace Grid; +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) + + LanczosParameters() { + ////////////////////////////// Default values + mass = 0; + ///////////////////////////////// + } + + 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(); + } + +}; + +} + +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? +} + + +typedef WilsonFermionD WilsonOp; +typedef typename WilsonFermionD::FermionField FermionField; + + +#if 0 // Hermitize a DWF operator by squaring it template class SquaredLinearOperator : public LinearOperatorBase { @@ -286,6 +351,7 @@ public: std::cout< void testSchurFromHess(Arnoldi& Arn, Field& src, int Nlarge, int Nm, int Nk) { @@ -332,21 +398,24 @@ int main (int argc, char ** argv) const int Ls=16; // GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); - std::vector lat_size {16, 16, 16, 32}; - std::cout << "Lattice size: " << lat_size << std::endl; - GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, +// std::vector lat_size {32, 32, 32, 32}; +// std::cout << "Lattice size: " << lat_size << std::endl; + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); +// GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); +// GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridCartesian * FGrid = UGrid; + GridRedBlackCartesian * FrbGrid = UrbGrid; // Construct a coarsened grid // poare TODO: replace this with the following line? - Coordinate clatt = lat_size; + Coordinate clatt = GridDefaultLatt(); // Coordinate clatt = GridDefaultLatt(); // [PO] initial line before I edited it for(int d=0;d PM; PM(PVdagM, src); + int Nm = 50; + int Nk = 12; + int Np = 38; + // int Nk = Nm+1; // if just running once + int maxIter = 10000; + int Nstop = 10; + RealD resid = 1.0e-5; + + std::vector boundary = {1,1,1,-1}; + WilsonOp::ImplParams Params(boundary); + +// DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); +// DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5); + + mass=LanParams.mass; + std::cout << GridLogIRL<< "mass "< PVdagM_t; - typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef ShiftedComplexPVdagMLinearOperator ShiftedComplexPVdagM_t; - PVdagM_t PVdagM(Ddwf, Dpv); - ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); - SquaredLinearOperator Dsq (Ddwf); - NonHermitianLinearOperator DLinOp (Ddwf); +// typedef PVdagMLinearOperator PVdagM_t; +// typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; +// typedef ShiftedComplexPVdagMLinearOperator ShiftedComplexPVdagM_t; +// PVdagM_t PVdagM(Ddwf, Dpv); +// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); +// SquaredLinearOperator Dsq (Ddwf); +// NonHermitianLinearOperator DLinOp (Ddwf); + + + NonHermitianLinearOperator Dwilson(WilsonOperator); /// <----- + MdagMLinearOperator HermOp(WilsonOperator); /// <----- + Gamma5HermitianLinearOperator HermOp2(WilsonOperator); /// <---- // PowerMethod PM; PM(PVdagM, src); - int Nm = 10; - int Nk = 4; - // int Nk = Nm+1; // if just running once - // int maxIter = 5; - // int maxIter = 1; - int maxIter = 5; - // int maxIter = 100; - int Nstop = 4; + + resid=LanParams.resid; + Nstop=LanParams.Nstop; + Nk=LanParams.Nk; + Np=LanParams.Np; + Nm = Nk + Np; + int Nu=16; + std::vector src(Nu,FGrid); + for(int i=0;i using namespace std; using namespace Grid; +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) + + LanczosParameters() { + ///////////////////////////////// + } + + template + LanczosParameters(Reader & TheReader){ + initialize(TheReader); + } + + template < class ReaderClass > + void initialize(Reader &TheReader){ + 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(); + } + +}; + +} + +#if 0 // Hermitize a DWF operator by squaring it template class SquaredLinearOperator : public LinearOperatorBase { @@ -288,6 +332,7 @@ public: std::cout< Date: Wed, 26 Nov 2025 22:13:27 +0000 Subject: [PATCH 2/2] Starting to modified KS --- Grid/algorithms/iterative/KrylovSchur.h | 34 +++++++++++++++++-------- examples/Example_krylov_schur.cc | 26 ++++++++++--------- 2 files changed, 37 insertions(+), 23 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index d5bf6a20..a83983dc 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -89,6 +89,14 @@ struct ComplexComparator RitzFilter RF; ComplexComparator (RitzFilter _rf) : RF(_rf) {} bool operator()(std::complex z1, std::complex z2) { + RealD tmp1, tmp2; + tmp1=std::abs(std::imag(z1)); + tmp2=std::abs(std::imag(z2)); +// if ( std::abs(std::real(z1)) <0.76) tmp1 +=4.; +// if ( std::abs(std::real(z2)) <0.76) tmp2 +=4.; + if ( std::abs(std::real(z1)) >2.) tmp1 +=4.; + if ( std::abs(std::real(z2)) >2.) tmp2 +=4.; + switch (RF) { case EvalNormSmall: return std::abs(z1) < std::abs(z2); @@ -103,9 +111,11 @@ struct ComplexComparator case EvalImLarge: return std::imag(z1) > std::imag(z2); case EvalImNormSmall: - return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + return tmp1 < tmp2; +// return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); case EvalImNormLarge: - return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); + return tmp1 > tmp2; +// return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); default: assert(0); } @@ -224,8 +234,8 @@ class ComplexSchurDecomposition { */ // void schurReorder(int Nk, std::function compare) { int schurReorder(int Nk) { - for (int i = 0; i < Nm; i++) { - for (int k = 0; k <= Nm - 2; k++) { + for (int i = 0; i < Nk; i++) { + for (int k = 0; k <= Nm - 2-i ; k++) { int idx = Nm - 2 - k; // TODO use RitzFilter enum here // if (std::abs(S(idx, idx)) < std::abs(S(idx+1, idx+1))) { // sort by largest modulus of eigenvalue @@ -235,10 +245,10 @@ class ComplexSchurDecomposition { } } } - int temp=1; - for (int i = 0; i < Nm-1; i++) { + int temp=Nk/2; + for (int i = temp; i < Nk; i++) { std::cout << GridLogMessage << "S " << i << " "<= Nstop || i == MaxIter - 1) { std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " @@ -469,7 +482,6 @@ class KrylovSchur { u = Zero(); } else { std::cout << GridLogMessage << "start= " < src(Nu,FGrid); for(int i=0;i