Merge branch 'feature/distil' of github.com:mmphys/Grid into feature/distil

Hadrons/Modules/MDistil/PerambFromSolve.hpp

@@ -42,6 +42,8 @@ public:
                                     std::string, PerambFileName,
                                     std::string, solve,
                                     int, nvec,
+                                    int, nvec_reduced,
+                                    int, LI_reduced,
                                     DistilParameters, Distil);
 };
 
@@ -125,8 +127,12 @@ void TPerambFromSolve<FImpl>::setup(void)
         grid3d = MakeLowerDimGrid(grid4d);
 
 
+	const int nvec_reduced{par().nvec_reduced};
+	const int LI_reduced{par().LI_reduced};
+    //envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
+    //          sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
     envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
-              sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
+              sIndexNames,Distil.Nt,nvec_reduced,LI_reduced,Distil.nnoise,Distil.Nt_inv,Distil.SI);
     envCreate(std::vector<Complex>, getName() + "_noise", 1,
               nvec*Distil.Ns*Distil.Nt*Distil.nnoise);
 
@@ -152,6 +158,8 @@ void TPerambFromSolve<FImpl>::Cleanup(void)
 template <typename FImpl>
 void TPerambFromSolve<FImpl>::execute(void)
 {
+	const int nvec_reduced{par().nvec_reduced};
+	const int LI_reduced{par().LI_reduced};
     const int nvec{par().nvec};
     const DistilParameters & Distil{par().Distil};
     const int LI{Distil.LI};
@@ -181,14 +189,14 @@ void TPerambFromSolve<FImpl>::execute(void)
 
 
     for (int inoise = 0; inoise < nnoise; inoise++) {
-      for (int dk = 0; dk < LI; dk++) {
+      for (int dk = 0; dk < LI_reduced; dk++) {
         for (int dt = 0; dt < Nt_inv; dt++) {
           for (int ds = 0; ds < Ns; ds++) {
             for (int is = 0; is < Ns; is++) {
               result_nospin = peekSpin(solve[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))],is);
               for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) {
                 ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Grid::QCD::Tdir);
-                for (int ivec = 0; ivec < nvec; ivec++) {
+                for (int ivec = 0; ivec < nvec_reduced; ivec++) {
                   ExtractSliceLocal(evec3d,epack.evec[ivec],0,t,3);
                   pokeSpin(perambulator(t, ivec, dk, inoise,dt,ds),innerProduct(evec3d, result_3d),is);
                   std::cout <<  "perambulator(t, ivec, dk, inoise,dt,ds)(is) = (" << t << "," << ivec << "," << dk << "," << inoise << "," << dt << "," << ds << ")(" << is << ") = " <<  perambulator(t, ivec, dk, inoise,dt,ds)()(is)() << std::endl;

Hadrons/Modules/MDistil/PerambMultipleSolves.hpp

@@ -25,9 +25,6 @@ class PerambMultipleSolvesPar: Serializable
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(PerambMultipleSolvesPar,
 		                    std::string, eigenPack,
-                                    std::string, PerambFileName,
-                                    std::string, ConfigFileDir,
-                                    std::string, ConfigFileName,
                                     std::string, UniqueIdentifier,
 				    int, nsolves,
 				    std::vector<int>, nvecs,
@@ -102,6 +99,7 @@ std::vector<std::string> TPerambMultipleSolves<FImpl>::getOutput(void)
       out.push_back(getName()+ "_solve_" +std::to_string(nvecs[i]));
     }
     
+      out.push_back(getName()+ "_noise");
     return out;
 }
 
@@ -112,7 +110,6 @@ void TPerambMultipleSolves<FImpl>::setup(void)
     Cleanup();
 
     const int nvec{par().nvec};
-    // auto &noise = envGet(std::vector<std::vector<std::vector<SpinVector>>>, par().noise);
     const int nsolves{par().nsolves};
     std::vector<int> nvecs{par().nvecs};
     const DistilParameters & Distil{par().Distil};
@@ -133,14 +130,11 @@ void TPerambMultipleSolves<FImpl>::setup(void)
 
     envTmpLat(GaugeField, "Umu");
     envTmpLat(LatticeSpinColourVector, "dist_source");
-    envTmp(std::vector<LatticeSpinColourVector>, "sources", 1,
+    //envTmp(std::vector<LatticeSpinColourVector>, "sources", 1,
-           std::vector<LatticeSpinColourVector>( nsolves, grid4d ));
+    //       std::vector<LatticeSpinColourVector>( nsolves, grid4d ));
     envTmpLat(LatticeSpinColourVector, "tmp2");
     envTmpLat(LatticeSpinColourVector, "result");
-    //envTmpLat(LatticeSpinColourVector, "result_single_component");
     envTmpLat(LatticeColourVector, "result_nospin");
-    //envTmpLat(LatticeColourVector, "tmp_nospin");
-    //envTmpLat(LatticeSpinVector, "peramb_tmp");
     envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
     envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
     envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
@@ -171,7 +165,6 @@ void TPerambMultipleSolves<FImpl>::execute(void)
     const int nvec{par().nvec};
     std::vector<int> nvecs{par().nvecs};
     const DistilParameters & Distil{par().Distil};
-    //const SolverParameters & Solver{par().Solver};
     const int LI{Distil.LI};
     //const int SI{Distil.SI};
     const int TI{Distil.TI};
@@ -190,8 +183,6 @@ void TPerambMultipleSolves<FImpl>::execute(void)
     const bool full_tdil{TI==Nt};
     const bool exact_distillation{full_tdil && LI==nvec};
 
-    const std::string &ConfigFileDir{par().ConfigFileDir};
-    const std::string &ConfigFileName{par().ConfigFileName};
     const std::string &UniqueIdentifier{par().UniqueIdentifier};
 
     auto        &noise   = envGet(std::vector<Complex>, getName() + "_noise");
@@ -203,44 +194,8 @@ void TPerambMultipleSolves<FImpl>::execute(void)
       solves[i]=unsmeared_sink;
     }
 
-    envGetTmp(GaugeField, Umu);
-  FieldMetaData header;
-  if((1)){
-    const std::vector<int> seeds({1, 2, 3, 4, 5});
-    GridParallelRNG pRNG4d(grid4d);
-    pRNG4d.SeedFixedIntegers(seeds);
-    std::cout << GridLogMessage << "now hot config" << std::endl;
-    SU<Nc>::HotConfiguration(pRNG4d, Umu);
-    std::cout << GridLogMessage << "hot cfg done." << std::endl;
-
-  // Set up the SAME gauge field on every time plane
-  //  int Nt = grid4d->gDimensions()[Tdir];
-  Grid_unquiesce_nodes();
-  
-  auto Usft = Umu;
-  Lattice<iScalar<vInteger> > coor(grid4d);
-  LatticeCoordinate(coor,Tdir);
-  for(int t=1;t<Nt;t++){
-    // t=1
-    //  Umu                Usft
-    // 0,1,2,3,4,5,6,7 -> 7,0,1,2,3,4,5,6 t=1
-    // 0,0,2,3,4,5,6,7    6,7,0,1,2,3,4,5 t=2
-    // 0,0,0,3,4,5,6,7    5,6,7,0,1,2,3,4 t=3
-    //...
-    
-    Usft = Cshift(Usft,Tdir,-1);
-    Umu = where(coor==t,Usft,Umu);
-  }
-  } else {
-    //std::string   fileName( "/home/dp008/dp008/dc-rich6/Scripts/ConfigsDeflQED/ckpoint_lat.3000" );
-    std::string   fileName(ConfigFileDir + ConfigFileName);
-    std::cout << GridLogMessage << "Loading NERSC configuration from '" << fileName << "'" << std::endl;
-    NerscIO::readConfiguration(Umu, header, fileName);
-    std::cout << GridLogMessage << "reading done." << std::endl;
-  }
-
     GridSerialRNG sRNG; 
-    sRNG.SeedUniqueString(ConfigFileName + "_"  + UniqueIdentifier);
+    sRNG.SeedUniqueString(UniqueIdentifier + std::to_string(vm().getTrajectory()));
     Real rn;
     
     for (int inoise=0;inoise<nnoise;inoise++) {
@@ -269,15 +224,11 @@ void TPerambMultipleSolves<FImpl>::execute(void)
   envGetTmp(LatticeColourVector, tmp3d_nospin);
   envGetTmp(LatticeColourVector, result_3d);
   envGetTmp(LatticeColourVector, evec3d);
-  envGetTmp(std::vector<LatticeSpinColourVector>, sources);
+  //envGetTmp(std::vector<LatticeSpinColourVector>, sources);
 
     const int Ntlocal{grid4d->LocalDimensions()[3]};
     const int Ntfirst{grid4d->LocalStarts()[3]};
 
-    std::cout << "init RBG "  << std::endl;
-    GridRedBlackCartesian RBGrid(grid4d);
-    std::cout << "init RBG done"  << std::endl;
-   
     {
 
     int t_inv;
@@ -287,10 +238,6 @@ void TPerambMultipleSolves<FImpl>::execute(void)
           for (int ds = 0; ds < Ns; ds++) {
             std::cout <<  "LapH source vector from noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl;
             dist_source = zero;
-            for (int isource = 0; isource < nsolves; isource ++){ 
-              if(dk < nvecs[isource])
-	        sources[isource] = zero;
-	    }
             tmp3d_nospin = zero;
             evec3d = zero;
             for (int it = dt; it < Nt; it += TI){
@@ -305,19 +252,13 @@ void TPerambMultipleSolves<FImpl>::execute(void)
                     tmp2=zero;
                     InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Grid::QCD::Tdir);
                     dist_source += tmp2;
-                    for (int isource = 0; isource < nsolves; isource ++){ 
-                      if(dk < nvecs[isource])
-		        sources[isource] += tmp2;
-		    }
                   }
                 }
               }
             }
             std::cout <<  "Inversion for noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl;
             result=zero;
-            for (int isource = 0; isource < nsolves; isource ++){ 
+	    v4dtmp = dist_source;
-              if(dk < nvecs[isource]){
-	        v4dtmp = sources[isource];
 	    if (Ls_ == 1){
 	      solver(result, v4dtmp);
     	    } else {
@@ -326,6 +267,8 @@ void TPerambMultipleSolves<FImpl>::execute(void)
 	       mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp);
 	       result = v4dtmp;
 	    }
+            for (int isource = 0; isource < nsolves; isource ++){ 
+              if(dk < nvecs[isource]){
                 solves[isource][inoise+nnoise*(dk+nvecs[isource]*(dt+Nt_inv*ds))] = result;
 	      }
 	    }
@@ -335,6 +278,9 @@ void TPerambMultipleSolves<FImpl>::execute(void)
     }
   }
 }
+
+
+
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE

tests/hadrons/Test_hadrons_distil.cc

@@ -87,10 +87,6 @@ void test_LapEvec(Application &application)
   application.createModule<MGauge::Random>(szGaugeName);
   // Now make an instance of the LapEvec object
   MDistil::LapEvecPar p;
-  //p.ConfigFileDir="/home/dp008/dp008/dc-rich6/Scripts/ConfigsDeflQED/";
-  //p.ConfigFileName="ckpoint_lat.3000";
-  p.ConfigFileDir="/home/dp008/dp008/paboyle/A2A/run/";
-  p.ConfigFileName="ckpoint_lat.IEEE64BIG.1100";
   p.gauge = szGaugeName;
   //p.EigenPackName = "ePack";
   //p.Distil.TI = 8;
@@ -120,8 +116,6 @@ void test_Perambulators(Application &application)
   MDistil::PerambLight::Par PerambPar;
   PerambPar.eigenPack="LapEvec";
   PerambPar.PerambFileName="peramb.bin";
-  PerambPar.ConfigFileDir="/home/dp008/dp008/paboyle/A2A/run/";
-  PerambPar.ConfigFileName="ckpoint_lat.IEEE64BIG.1100";
   PerambPar.UniqueIdentifier="full_dilution";
   PerambPar.solver="CG_s";
   PerambPar.Distil.tsrc = 0;
@@ -141,6 +135,99 @@ void test_Perambulators(Application &application)
   application.createModule<MDistil::PerambLight>("Peramb",PerambPar);
 }
 /////////////////////////////////////////////////////////////
+// Multiple Perambulators
+/////////////////////////////////////////////////////////////
+
+void test_MultiPerambulators(Application &application)
+{
+  // PerambLight parameters
+  MDistil::PerambLight::Par PerambPar;
+  PerambPar.eigenPack="LapEvec";
+  PerambPar.UniqueIdentifier="full_dilution";
+  PerambPar.PerambFileName="Peramb5";
+  PerambPar.solver="CG_s";
+  PerambPar.Distil.tsrc = 0;
+  PerambPar.Distil.nnoise = 1;
+  PerambPar.Distil.LI=5;
+  PerambPar.Distil.SI=4;
+  PerambPar.Distil.TI=8;
+  PerambPar.nvec=5;
+  PerambPar.Distil.Ns=4;
+  PerambPar.Distil.Nt=8;
+  PerambPar.Distil.Nt_inv=1;
+  application.createModule<MDistil::PerambLight>("Peramb5",PerambPar);
+  MDistil::PerambFromSolve::Par SolvePar;
+  SolvePar.eigenPack="LapEvec";
+  SolvePar.PerambFileName="Peramb2";
+  SolvePar.solve = "Peramb5_unsmeared_sink";
+  SolvePar.Distil.nnoise = 1;
+  SolvePar.Distil.LI=5;
+  SolvePar.Distil.SI=4;
+  SolvePar.Distil.TI=8;
+  SolvePar.nvec=5;
+  SolvePar.nvec_reduced=2;
+  SolvePar.LI_reduced=2;
+  SolvePar.Distil.Ns=4;
+  SolvePar.Distil.Nt=8;
+  SolvePar.Distil.Nt_inv=1;
+  application.createModule<MDistil::PerambFromSolve>("Peramb2",SolvePar);
+  SolvePar.PerambFileName="Peramb3";
+  SolvePar.nvec_reduced=3;
+  SolvePar.LI_reduced=3;
+  application.createModule<MDistil::PerambFromSolve>("Peramb3",SolvePar);
+  MDistil::DistilVectors::Par DistilVecPar;
+  DistilVecPar.noise="Peramb5_noise";
+  DistilVecPar.perambulator="Peramb2";
+  DistilVecPar.eigenPack="LapEvec";
+  DistilVecPar.tsrc = 0;
+  DistilVecPar.nnoise = 1;
+  DistilVecPar.LI=2;
+  DistilVecPar.SI=4;
+  DistilVecPar.TI=8;
+  DistilVecPar.nvec=2;
+  DistilVecPar.Ns=4;
+  DistilVecPar.Nt=8;
+  DistilVecPar.Nt_inv=1;
+  application.createModule<MDistil::DistilVectors>("DistilVecs2",DistilVecPar);
+  DistilVecPar.perambulator="Peramb3";
+  DistilVecPar.LI=3;
+  DistilVecPar.nvec=3;
+  application.createModule<MDistil::DistilVectors>("DistilVecs3",DistilVecPar);
+  DistilVecPar.perambulator="Peramb5_perambulator_light";
+  DistilVecPar.LI=5;
+  DistilVecPar.nvec=5;
+  application.createModule<MDistil::DistilVectors>("DistilVecs5",DistilVecPar);
+  MContraction::A2AMesonField::Par A2AMesonFieldPar;
+  A2AMesonFieldPar.left="DistilVecs2_rho";
+  A2AMesonFieldPar.right="DistilVecs2_rho";
+  A2AMesonFieldPar.output="MesonSinksRho2";
+  A2AMesonFieldPar.gammas="all";
+  A2AMesonFieldPar.mom={"0 0 0"};
+  A2AMesonFieldPar.cacheBlock=2;
+  A2AMesonFieldPar.block=4;
+  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho2",A2AMesonFieldPar);
+  A2AMesonFieldPar.left="DistilVecs2_phi";
+  A2AMesonFieldPar.right="DistilVecs2_phi";
+  A2AMesonFieldPar.output="MesonSinksPhi2";
+  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi2",A2AMesonFieldPar);
+  A2AMesonFieldPar.left="DistilVecs3_rho";
+  A2AMesonFieldPar.right="DistilVecs3_rho";
+  A2AMesonFieldPar.output="MesonSinksRho3";
+  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho3",A2AMesonFieldPar);
+  A2AMesonFieldPar.left="DistilVecs3_phi";
+  A2AMesonFieldPar.right="DistilVecs3_phi";
+  A2AMesonFieldPar.output="MesonSinksPhi3";
+  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi3",A2AMesonFieldPar);
+  A2AMesonFieldPar.left="DistilVecs5_rho";
+  A2AMesonFieldPar.right="DistilVecs5_rho";
+  A2AMesonFieldPar.output="MesonSinksRho5";
+  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho5",A2AMesonFieldPar);
+  A2AMesonFieldPar.left="DistilVecs5_phi";
+  A2AMesonFieldPar.right="DistilVecs5_phi";
+  A2AMesonFieldPar.output="MesonSinksPhi5";
+  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi5",A2AMesonFieldPar);
+}
+/////////////////////////////////////////////////////////////
 // DistilVectors
 /////////////////////////////////////////////////////////////
 
@@ -168,8 +255,6 @@ void test_PerambulatorsS(Application &application)
   MDistil::PerambLight::Par PerambPar;
   PerambPar.eigenPack="LapEvec";
   PerambPar.PerambFileName="perambS.bin";
-  PerambPar.ConfigFileDir="/home/dp008/dp008/paboyle/A2A/run/";
-  PerambPar.ConfigFileName="ckpoint_lat.IEEE64BIG.1100";
   PerambPar.UniqueIdentifier="full_dilution";
   PerambPar.solver="CG_s";
   PerambPar.Distil.tsrc = 0;
@@ -991,6 +1076,12 @@ int main(int argc, char *argv[])
       test_DistilVectorsAslashSeq( application );
       test_MesonFieldAslashSeq( application );
       break;
+    case 13:
+      test_Global( application );
+      test_SolverS( application );
+      test_LapEvec( application );
+      test_MultiPerambulators( application );
+      break;
   }
   LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;