Checking in before pulling

Grid/algorithms/iterative/Arnoldi.h

@@ -51,14 +51,23 @@ struct ComplexComparator
   ComplexComparator (RitzFilter _f) : f(_f) {}
   bool operator()(std::complex<double> z1, std::complex<double> 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);
     }

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) {
+    int schurReorder(int Nk) {
-      for (int i = 0; i < Nk; i++) {
+      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 << " "<<std::real(S(i,i))<<" "<<std::imag(S(i,i)) << std::endl;
+	if (std::abs(std::imag(S(i,i))< 0.01)) temp= i+1;
+      }
+      return temp;
     }
 
     void schurReorderBlock() {
@@ -197,7 +202,6 @@ class KrylovSchur {
 
     Eigen::VectorXcd evals;                 // evals of Rayleigh quotient
     Eigen::MatrixXcd littleEvecs;           // Nm x Nm evecs matrix
-    std::vector<Field> evecs;               // Vector of evec fields
 
     RitzFilter ritzFilter;                        // how to sort evals
 
@@ -210,6 +214,8 @@ class KrylovSchur {
       u = Zero();
     };
 
+    std::vector<Field> 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= "<<Nk<<std::endl;
         Eigen::MatrixXcd Q = schur.getMatrixQ();
         Qt = Q.adjoint();                           // TODO should Q be real?
         Eigen::MatrixXcd S = schur.getMatrixS();
@@ -286,7 +294,7 @@ class KrylovSchur {
 
         // checkKSDecomposition();
 
-        std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** " << std::endl;
+        std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** Nk=" << Nk << std::endl;
 
         std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl;
 
@@ -354,7 +362,8 @@ class KrylovSchur {
         Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm);
         u = Zero();
       } else {
-        assert( start == basis.size() );      // should be starting at the end of basis (start = Nk)
+        std::cout << GridLogMessage << "start= " <<start<< " basis.size= "<<basis.size()<<std::endl;
+        assert( start <= basis.size() );      // should be starting at the end of basis (start = Nk)
         std::cout << GridLogMessage << "Resetting Rayleigh and b" << std::endl;
         Eigen::MatrixXcd RayleighCp = Rayleigh;
         Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm);
@@ -376,7 +385,7 @@ class KrylovSchur {
         }
 
         if (doubleOrthog) {
-          std::cout << GridLogMessage << "Double orthogonalizing." << std::endl;
+//          std::cout << GridLogMessage << "Double orthogonalizing." << std::endl;
           for (int j = 0; j < basis.size(); j++) {
             coeff = innerProduct(basis[j], w);      // see if there is any residual component in basis[j] direction
             Rayleigh(j, i) += coeff;                // if coeff is non-zero, adjust Rayleigh

examples/Example_krylov_schur.cc

@@ -41,6 +41,71 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 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 <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();
+  }
+  
+};
+
+}
+
+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?
+}
+
+
+typedef WilsonFermionD WilsonOp;
+typedef typename WilsonFermionD::FermionField FermionField;
+
+
+#if 0
 // Hermitize a DWF operator by squaring it
 template<class Matrix,class Field>
 class SquaredLinearOperator : public LinearOperatorBase<Field> {
@@ -286,6 +351,7 @@ public:
     std::cout<<GridLogMessage << "Done " <<std::endl;
   }
 };
+#endif
 
 template<class Field>
 void testSchurFromHess(Arnoldi<Field>& 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<int> lat_size {16, 16, 16, 32};
+//  std::vector<int> lat_size {32, 32, 32, 32};
-  std::cout << "Lattice size: " << lat_size << std::endl;
+//  std::cout << "Lattice size: " << lat_size << std::endl;
-  GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, 
+  GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
 								          GridDefaultSimd(Nd,vComplex::Nsimd()),
 								          GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
 
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+//  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(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<clatt.size();d++){
+  std::cout << GridLogMessage<< clatt[d] <<std::endl;
     clatt[d] = clatt[d]/2;
     //    clatt[d] = clatt[d]/4;
   }
@@ -360,7 +429,6 @@ int main (int argc, char ** argv)
   GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
   GridParallelRNG          CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
 
-  LatticeFermion    src(FGrid); random(RNG5,src);
   LatticeFermion result(FGrid); result=Zero();
   LatticeFermion    ref(FGrid); ref=Zero();
   LatticeFermion    tmp(FGrid);
@@ -369,14 +437,56 @@ int main (int argc, char ** argv)
 
   FieldMetaData header;
 //   std::string file("ckpoint_lat.4000");
-  std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_lat.4000");
+//  std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_lat.4000");
+  std::string file("./config");
   NerscIO::readConfiguration(Umu,header,file);
 
+  LanczosParameters LanParams;
+  {
+    XmlReader  HMCrd("LanParams.xml");
+    read(HMCrd,"LanczosParameters",LanParams);
+  }
+
+  std::cout << GridLogMessage<< LanParams <<std::endl;
+  {
+    XmlWriter HMCwr("LanParams.xml.out");
+    write(HMCwr,"LanczosParameters",LanParams);
+  }
+
+#if 0
+  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();
+  }
+#endif
+
+
   RealD mass=0.01;
   RealD M5=1.8;
 
-  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  // PowerMethod<LatticeFermion> PM; PM(PVdagM, src);
-  DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
+  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<Complex> 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 "<<mass<<std::endl;
+  WilsonOp WilsonOperator(Umu,*UGrid,*UrbGrid,mass,Params);
 
   // const int nbasis = 20;            // size of approximate basis for low-mode space
   const int nbasis = 3;            // size of approximate basis for low-mode space
@@ -392,34 +502,61 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
   std::cout<<GridLogMessage<<std::endl;
 
-  typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
+//  typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
-  typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
+//  typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
-  typedef ShiftedComplexPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedComplexPVdagM_t;
+//  typedef ShiftedComplexPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedComplexPVdagM_t;
-  PVdagM_t PVdagM(Ddwf, Dpv);
+//  PVdagM_t PVdagM(Ddwf, Dpv);
-  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv);
+//  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv);
-  SquaredLinearOperator<DomainWallFermionD, LatticeFermionD> Dsq (Ddwf);
+//  SquaredLinearOperator<DomainWallFermionD, LatticeFermionD> Dsq (Ddwf);
-  NonHermitianLinearOperator<DomainWallFermionD, LatticeFermionD> DLinOp (Ddwf);
+//  NonHermitianLinearOperator<DomainWallFermionD, LatticeFermionD> DLinOp (Ddwf);
+
+
+  NonHermitianLinearOperator<WilsonOp,FermionField> Dwilson(WilsonOperator); /// <-----
+  MdagMLinearOperator<WilsonOp,FermionField> HermOp(WilsonOperator); /// <-----
+  Gamma5HermitianLinearOperator <WilsonOp,LatticeFermion> HermOp2(WilsonOperator); /// <----
 
   // PowerMethod<LatticeFermion> PM; PM(PVdagM, src);
-  int Nm = 10;
+
-  int Nk = 4; 
+  resid=LanParams.resid;
-  // int Nk = Nm+1;     // if just running once
+  Nstop=LanParams.Nstop;
-  // int maxIter = 5;
+  Nk=LanParams.Nk;
-  // int maxIter = 1;
+  Np=LanParams.Np;
-  int maxIter = 5;
+  Nm = Nk + Np;
-  // int maxIter = 100;
+  int Nu=16;
-  int Nstop = 4;
+  std::vector<LatticeFermion> src(Nu,FGrid); 
+  for(int i=0;i<Nu;i++) random(RNG5,src[i]);
 
   Coordinate origin ({0,0,0,0});
-  auto tmpSrc = peekSite(src, origin);
+  auto tmpSrc = peekSite(src[0], origin);
   std::cout << "[DEBUG] Source at origin = " <<  tmpSrc << std::endl;
-  LatticeFermion src2 = src;
+  LatticeFermion src2 = src[0];
 
   // Run KrylovSchur and Arnoldi on a Hermitian matrix
-  std::cout << GridLogMessage << "Runnning Krylov Schur" << std::endl;
+  std::cout << GridLogMessage << "Running Krylov Schur" << std::endl;
   // KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge);
-  KrylovSchur KrySchur (Dsq, FGrid, 1e-8);
+//  KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall);
-  KrySchur(src, maxIter, Nm, Nk, Nstop);
+//  Hacked, really EvalImagSmall
+  KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalReSmall);
+//  BlockKrylovSchur KrySchur (HermOp2, UGrid, Nu, resid,EvalNormSmall);
+  KrySchur(src[0], maxIter, Nm, Nk, Nstop);
+  std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl;
+
+  src[0]=KrySchur.evecs[0];
+  for (int i=1;i<Nstop;i++) src[0]+=KrySchur.evecs[i];
+  for (int i=0;i<Nstop;i++) 
+  {
+	std::string evfile ("./evec_"+std::to_string(mass)+"_"+std::to_string(i));
+        auto evdensity = localInnerProduct(KrySchur.evecs[i],KrySchur.evecs[i] );
+        writeFile(evdensity,evfile);
+
+  }
+
+  {
+        std::string evfile ("./evec_"+std::to_string(mass)+"_sum");
+//        auto evdensity = localInnerProduct(evec[i],evec[i] );
+        writeFile(src[0],evfile);
+  }
+
 
   /*
   std::cout << GridLogMessage << "Running Arnoldi" << std::endl;

examples/Example_spec_kryschur.cc

@@ -43,6 +43,50 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 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 <class ReaderClass >
+  LanczosParameters(Reader<ReaderClass> & TheReader){
+    initialize(TheReader);
+  }
+
+  template < class ReaderClass > 
+  void initialize(Reader<ReaderClass> &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 Matrix,class Field>
 class SquaredLinearOperator : public LinearOperatorBase<Field> {
@@ -288,6 +332,7 @@ public:
     std::cout<<GridLogMessage << "Done " <<std::endl;
   }
 };
+#endif
 
 int main (int argc, char ** argv)
 {