Checking in before pulling KrylovSchur

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<b*b) return true;
+	return false;
+}
+#endif
+
 template<class Field> class ImplicitlyRestartedLanczosHermOpTester  : public ImplicitlyRestartedLanczosTester<Field>
 {
  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<Field> & 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<RealD>());
       std::cout<<GridLogIRL <<" evals sorted "<<std::endl;
+//      eval2_copy = eval2;
       const int chunk=8;
       for(int io=0; io<k2;io+=chunk){
 	std::cout<<GridLogIRL << "eval "<< std::setw(3) << io ;
@@ -333,6 +352,7 @@ until convergence
       //////////////////////////////////
       Qt = Eigen::MatrixXd::Identity(Nm,Nm);
       for(int ip=k2; ip<Nm; ++ip){ 
+//        std::cout<<GridLogIRL <<"QR decompose "<<eval2[ip]<<std::endl;
 	QR_decomp(eval,lme,Nm,Nm,Qt,eval2[ip],k1,Nm);
       }
       std::cout<<GridLogIRL <<"QR decomposed "<<std::endl;
@@ -375,7 +395,8 @@ until convergence
 
 	//  power of two search pattern;  not every evalue in eval2 is assessed.
 	int allconv =1;
-	for(int jj = 1; jj<=Nstop; jj*=2){
+//	for(int jj = 1; jj<=Nstop; jj*=2){
+	for(int jj = 1; jj<=Nstop; jj++){
 	  int j = Nstop-jj;
 	  RealD e = eval2_copy[j]; // Discard the evalue
 	  basisRotateJ(B,evec,Qt,j,0,Nk,Nm);	    

tests/lanczos/Test_wilson_lanczos.cc

@@ -6,7 +6,7 @@ Source file: ./tests/Test_dwf_lanczos.cc
 
 Copyright (C) 2015
 
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+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
@@ -27,6 +27,9 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid/Grid.h>
+#include <Grid/parallelIO/IldgIOtypes.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
+
 
 using namespace std;
 using namespace Grid;
@@ -38,18 +41,111 @@ typedef typename WilsonFermionD::FermionField FermionField;
 
 RealD AllZero(RealD x) { return 0.; }
 
+template <class T> 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 <class ReaderClass >
+  LanczosParameters(Reader<ReaderClass> & TheReader){
+    initialize(TheReader);
+  }
+
+  template < class ReaderClass > 
+  void initialize(Reader<ReaderClass> &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<int> nblock(4,1);
+  std::vector<int> mpi_split(4,1);
+//Interested in avoiding degeneracy only for now
+  nblock[3]=2;
+
+  int mrhs=1;
+  for(int i =0;i<Ndir;i++){
+      mpi_split[i] = mpi_layout[i] / nblock[i];
+      mrhs *= nblock[i];
+  }
+
+
   GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
       GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
       GridDefaultMpi());
-  GridRedBlackCartesian* UrbGrid =
+
-      SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian * SGrid = new GridCartesian(GridDefaultLatt(),
+                                                    GridDefaultSimd(Nd,vComplex::Nsimd()),
+                                                    mpi_split,
+                                                    *UGrid);
+
+  GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
   GridCartesian* FGrid = UGrid;
   GridRedBlackCartesian* FrbGrid = UrbGrid;
-  printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid,
+//  printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid, FrbGrid);
-         FrbGrid);
 
   std::vector<int> seeds4({1, 2, 3, 4});
   std::vector<int> seeds5({5, 6, 7, 8});
@@ -62,12 +158,15 @@ int main(int argc, char** argv) {
 
   LatticeGaugeField Umu(UGrid);
 //  SU<Nc>::HotConfiguration(RNG4, Umu);
-  SU<Nc>::ColdConfiguration(Umu);
+//  SU<Nc>::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<LatticeColourMatrix> U(4, UGrid);
@@ -75,38 +174,101 @@ int main(int argc, char** argv) {
     U[mu] = PeekIndex<LorentzIndex>(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 <<std::endl;
+  { 
+    XmlWriter HMCwr("LanParams.xml.out");
+    write(HMCwr,"LanczosParameters",LanParams);
+  }
+
+  mass=LanParams.mass;
+  resid=LanParams.resid;
+  Nstop=LanParams.Nstop;
+  Nu = mrhs;
+  Nk=LanParams.Nk;
+  Np=LanParams.Np;
+  Nm = Nk + Np;
+
+//  FermionField src(FGrid);
+  std::vector<FermionField> src(Nu,FGrid);
+  for(int i =0;i<Nu;i++) gaussian(RNG5, src[i]);
+
+  if(LanParams.ReadEvec) {
+    std::string evecs_file="evec_in";
+    std::cout << GridLogIRL<< "Reading evecs from "<<evecs_file<<std::endl;
+    emptyUserRecord record;
+    Grid::ScidacReader RD;
+    RD.open(evecs_file);
+    RD.readScidacFieldRecord(src[0],record);
+    RD.close();
+  }
+
   std::vector<Complex> boundary = {1,1,1,-1};
 //  std::vector<Complex> 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<FermionOp,LatticeFermion> HermOp(WilsonOperator); /// <-----
+  MdagMLinearOperator<FermionOp,FermionField> HermOp(WilsonOperator); /// <-----
+  FermionOp WilsonSplit(s_Umu,*SFGrid,*SFrbGrid,mass,Params);
+  MdagMLinearOperator<FermionOp,FermionField> SHermOp(WilsonSplit); /// <-----
   //SchurDiagTwoOperator<FermionOp,FermionField> HermOp(WilsonOperator);
-
+  Gamma5HermitianLinearOperator <FermionOp,LatticeFermion> HermOp2(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;
 
   std::vector<double> Coeffs{0, 1.};
   Polynomial<FermionField> PolyX(Coeffs);
-  Chebyshev<FermionField> Cheby(0.0, 10., 12);
+//  Chebyshev<FermionField> Cheby(0.5, 60., 31);
+//                                  RealD, ChebyLow,
+//                                RealD, ChebyHigh,
+//                                Integer, ChebyOrder)
+
+  Chebyshev<FermionField> Cheby(LanParams.ChebyLow,LanParams.ChebyHigh,LanParams.ChebyOrder);
 
   FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
      PlainHermOp<FermionField> Op     (HermOp);
+     PlainHermOp<FermionField> Op2     (HermOp2);
 
-//  ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby, Op, Nstop, Nk, Nm, resid, MaxIt);
+//  ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby, Op2, Nstop, Nk, Nm, resid, MaxIt);
-  SimpleLanczos<FermionField> IRL(Op,Nstop, Nk, Nm, resid, MaxIt);
+//  SimpleLanczos<FermionField> IRL(Op,Nstop, Nk, Nm, resid, MaxIt);
-
+    ImplicitlyRestartedBlockLanczos<FermionField> IRBL(HermOp, SHermOp,
+                                                     FrbGrid,SFrbGrid,mrhs,
+                                                     Cheby,
+                                                     Nstop, Nstop*2,
+                                                     Nu, Nk, Nm,
+                                                     resid, MaxIt,
+                                                     IRBLdiagonaliseWithEigen);
+  IRBL.split_test=1;
   std::vector<RealD> eval(Nm);
-  FermionField src(FGrid);
-  gaussian(RNG5, src);
   std::vector<FermionField> 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<Nu;i++) src[i]=evec[i];
+  for(int i=Nu;i<Nstop;i++) src[i%Nu] +=evec[i];
+//  src  = evec[0]+evec[1]+evec[2];
+//  src  += evec[3]+evec[4]+evec[5];
+//  src  += evec[6]+evec[7]+evec[8];
+  mass += LanParams.mstep;
+}
 
   Grid_finalize();
 }

tests/lanczos/Test_wilson_specflow.cc

@@ -62,6 +62,8 @@ struct LanczosParameters: Serializable {
 				Integer, Nstop,
                                 Integer, Nk,
                                 Integer, Np,
+                                Integer, ReadEvec,
+	  			RealD, resid,
 	  			RealD, ChebyLow,
 	  			RealD, ChebyHigh,
 	  			Integer, ChebyOrder)
@@ -139,8 +141,8 @@ int main(int argc, char** argv) {
   FieldMetaData header;
   std::string file("./config");
 
-  int precision32 = 0;
+//  int precision32 = 0;
-  int tworow      = 0;
+//  int tworow      = 0;
 //  NerscIO::writeConfiguration(Umu,file,tworow,precision32);
   NerscIO::readConfiguration(Umu,header,file);
 
@@ -178,6 +180,7 @@ int main(int argc, char** argv) {
   }
 
   mass=LanParams.mass;
+  resid=LanParams.resid;
   Nstop=LanParams.Nstop;
   Nk=LanParams.Nk;
   Np=LanParams.Np;
@@ -185,12 +188,22 @@ int main(int argc, char** argv) {
 
   FermionField src(FGrid);
   gaussian(RNG5, src);
+  if(LanParams.ReadEvec) {
+    std::string evecs_file="evec_in";
+    std::cout << GridLogIRL<< "Reading evecs from "<<evecs_file<<std::endl;
+    emptyUserRecord record;
+    Grid::ScidacReader RD;
+    RD.open(evecs_file);
+    RD.readScidacFieldRecord(src,record);
+    RD.close();
+  }
+
   std::vector<Complex> boundary = {1,1,1,-1};
 //  std::vector<Complex> 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<FermionOp,FermionField> HermOp(WilsonOperator); /// <-----
   //SchurDiagTwoOperator<FermionOp,FermionField> 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 "<<evfile<<std::endl;
+    emptyUserRecord record;
+    Grid::ScidacReader RD;
+    RD.open(evfile);
+    RD.readScidacFieldRecord(evdensity,record);
+    RD.close();
+  }
     }
+    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);
   }
   src  = evec[0]+evec[1]+evec[2];
   src  += evec[3]+evec[4]+evec[5];