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

Conflicts: Grid/qcd/utils/BaryonUtils.h Hadrons/Modules/MContraction/Baryon2.hpp
2025-11-04 05:54:32 +00:00 · 2019-08-05 14:19:43 +01:00
parent 29df60c0cb 8d97e2a02a
commit b875edceab
8 changed files with 231 additions and 117 deletions
--- a/Grid/communicator/SharedMemoryMPI.cc
+++ b/Grid/communicator/SharedMemoryMPI.cc
@@ -143,8 +143,9 @@ void GlobalSharedMemory::OptimalCommunicator(const std::vector<int> &processors,
  gethostname(name,namelen);
  int nscan = sscanf(name,"r%di%dn%d",&R,&I,&N) ;

-  if(nscan==3) OptimalCommunicatorHypercube(processors,optimal_comm);
-  else         OptimalCommunicatorSharedMemory(processors,optimal_comm);
+  //  if(nscan==3) OptimalCommunicatorHypercube(processors,optimal_comm);
+  //  else         OptimalCommunicatorSharedMemory(processors,optimal_comm);
+  OptimalCommunicatorSharedMemory(processors,optimal_comm);
 }
 void GlobalSharedMemory::OptimalCommunicatorHypercube(const std::vector<int> &processors,Grid_MPI_Comm & optimal_comm)
 {
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -41,11 +41,8 @@ public:
 					 const Gamma GammaB,
 					 ComplexField &baryon_corr);
 
-// static PropagatorField quarkContract13(const PropagatorField &q1,
-//			     const PropagatorField &q2);
- static void quarkContract13(const PropagatorField &q1,
-			     const PropagatorField &q2,
-                             PropagatorField &q_out);
+ static LatticeSpinColourMatrix quarkContract13(const PropagatorField &q1,
+			     const PropagatorField &q2);
 };


@@ -319,31 +316,23 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1,
 //QDP / CHROMA - style diquark construction
 // (q_out)^{c'c}_{alpha,beta} = epsilon^{abc} epsilon^{a'b'c'} (q1)^{aa'}_{rho alpha}^* (q2)^{bb'}_{rho beta}
 template<class FImpl>
-//typename FImpl::PropagatorField BaryonUtils<FImpl>::quarkContract13(const PropagatorField &q1,
-//					 const PropagatorField &q2)
-void BaryonUtils<FImpl>::quarkContract13(const PropagatorField &q1,
-					 const PropagatorField &q2,
-                                         PropagatorField &q_out)
+LatticeSpinColourMatrix BaryonUtils<FImpl>::quarkContract13(const PropagatorField &q1,
+					 const PropagatorField &q2)
 {
  GridBase *grid = q1._grid;

-
  std::vector<std::vector<int>> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
  std::vector<int> epsilon_sgn = {1,1,1,-1,-1,-1};

-  //PropagatorField q_out;//=zero;
+  // TODO: Felix, made a few changes to fix this as there were compiler errors. Please validate!
+  LatticeSpinColourMatrix q_out(grid);
+  // q_out = zero; TODO: Don't think you need this, as you'll set each site explicitly anyway

  parallel_for(int ss=0;ss<grid->oSites();ss++){
-
-    typedef typename ComplexField::vector_object vobj;
-
-    auto D1 = q1._odata[ss];
-    auto D2 = q2._odata[ss];
-    //auto D_out = q_out._odata[ss];
-    //D_out=zero;
- 
-    pobj D_out=zero;
-
+    const auto & D1    = q1._odata[ss];
+    const auto & D2    = q2._odata[ss];
+          auto & D_out = q_out._odata[ss];
+    D_out=zero;
    for (int ie_src=0; ie_src < 6 ; ie_src++){
      int a_src = epsilon[ie_src][0]; //a
      int b_src = epsilon[ie_src][1]; //b
@@ -359,14 +348,10 @@ void BaryonUtils<FImpl>::quarkContract13(const PropagatorField &q1,
        }}}
      }
    }
-
-    q_out._odata[ss]=D_out;
-
 } //end loop over lattice sites

-  //return q_out;

+  return q_out;
 }

-
 }}
--- a/Hadrons/Modules.hpp
+++ b/Hadrons/Modules.hpp
@@ -60,7 +60,6 @@
 #include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
 #include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
 #include <Hadrons/Modules/MContraction/Baryon.hpp>
-#include <Hadrons/Modules/MContraction/Baryon_old.hpp>
 #include <Hadrons/Modules/MContraction/Meson.hpp>
 #include <Hadrons/Modules/MContraction/A2ALoop.hpp>
 #include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
--- a/Hadrons/Modules/MContraction/Baryon.hpp
+++ b/Hadrons/Modules/MContraction/Baryon.hpp
@@ -116,6 +116,12 @@ template <typename FImpl1, typename FImpl2, typename FImpl3>
 void TBaryon<FImpl1, FImpl2, FImpl3>::setup(void)
 {
    envTmpLat(LatticeComplex, "c");
+    envTmpLat(LatticeComplex, "c1");
+    envTmpLat(LatticeComplex, "c2");
+    envTmpLat(LatticeComplex, "c3");
+    envTmpLat(LatticeComplex, "c4");
+    envTmpLat(LatticeComplex, "c5");
+    envTmpLat(LatticeComplex, "c6");
    envTmpLat(LatticeComplex, "diquark");
  // Translate the full string naming the desired gamma structure into the one we need to use
  const std::string gamma{ par().gamma };
@@ -172,6 +178,12 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
    auto       &q2 = envGet(PropagatorField2, par().q2);
    auto       &q3 = envGet(PropagatorField3, par().q3);
    envGetTmp(LatticeComplex, c);
+    envGetTmp(LatticeComplex, c1);
+    envGetTmp(LatticeComplex, c2);
+    envGetTmp(LatticeComplex, c3);
+    envGetTmp(LatticeComplex, c4);
+    envGetTmp(LatticeComplex, c5);
+    envGetTmp(LatticeComplex, c6);
    envGetTmp(LatticeComplex, diquark);
    Result     result;
    int nt = env().getDim(Tp);
@@ -179,19 +191,63 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
    const std::string gamma{ par().gamma };
    std::vector<TComplex> buf;
    
+    Result     result1;
+    Result     result2;
+    Result     result3;
+    Result     result4;
+    Result     result5;
+    Result     result6;
+    result1.corr.resize(nt);
+    result2.corr.resize(nt);
+    result3.corr.resize(nt);
+    result4.corr.resize(nt);
+    result5.corr.resize(nt);
+    result6.corr.resize(nt);
+    std::vector<TComplex> buf1;
+    std::vector<TComplex> buf2;
+    std::vector<TComplex> buf3;
+    std::vector<TComplex> buf4;
+    std::vector<TComplex> buf5;
+    std::vector<TComplex> buf6;
+
    const Gamma GammaA{ Gamma::Algebra::Identity };
    const Gamma GammaB{ al };

-    BaryonUtils<FIMPL>::ContractBaryons(q1,q2,q3,GammaA,GammaB,c);
+    //BaryonUtils<FIMPL>::ContractBaryons(q1,q2,q3,GammaA,GammaB,c);
+    BaryonUtils<FIMPL>::ContractBaryons_debug(q1,q2,q3,GammaA,GammaB,c1,c2,c3,c4,c5,c6,c);

    sliceSum(c,buf,Tp);
+    sliceSum(c1,buf1,Tp);
+    sliceSum(c2,buf2,Tp);
+    sliceSum(c3,buf3,Tp);
+    sliceSum(c4,buf4,Tp);
+    sliceSum(c5,buf5,Tp);
+    sliceSum(c6,buf6,Tp);

    for (unsigned int t = 0; t < buf.size(); ++t)
    {
        result.corr[t] = TensorRemove(buf[t]);
+        result1.corr[t] = TensorRemove(buf1[t]);
+        result2.corr[t] = TensorRemove(buf2[t]);
+        result3.corr[t] = TensorRemove(buf3[t]);
+        result4.corr[t] = TensorRemove(buf4[t]);
+        result5.corr[t] = TensorRemove(buf5[t]);
+        result6.corr[t] = TensorRemove(buf6[t]);
    }

+    std::string ostr1{ par().output + "_1"};
+    std::string ostr2{ par().output + "_2"};
+    std::string ostr3{ par().output + "_3"};
+    std::string ostr4{ par().output + "_4"};
+    std::string ostr5{ par().output + "_5"};
+    std::string ostr6{ par().output + "_6"};
    saveResult(par().output, "baryon", result);
+    saveResult(ostr1, "baryon1", result1);
+    saveResult(ostr2, "baryon2", result2);
+    saveResult(ostr3, "baryon3", result3);
+    saveResult(ostr4, "baryon4", result4);
+    saveResult(ostr5, "baryon5", result5);
+    saveResult(ostr6, "baryon6", result6);
 }

 END_MODULE_NAMESPACE
--- a/Hadrons/Modules/MContraction/Baryon2.hpp
+++ b/Hadrons/Modules/MContraction/Baryon2.hpp
@@ -182,11 +182,13 @@ void TBaryon2<FImpl1, FImpl2, FImpl3>::execute(void)
    const Gamma GammaA{ Gamma::Algebra::Identity };
    const Gamma GammaB{ al };

-    //diquark = BaryonUtils<FIMPL>::quarkContract13(q2*GammaB,GammaB*q3);
-    BaryonUtils<FIMPL>::quarkContract13(q2*GammaB,GammaB*q3,diquark);
+    LatticeSpinColourMatrix diquark( q1._grid ); // TODO: Felix, I added "q1._grid". I presume this is correct?
+
+    diquark = BaryonUtils<FIMPL>::quarkContract13(q2*GammaB,GammaB*q3);

    //result = trace(GammaA*GammaA * traceColour(q1*traceSpin(diquark))) + 2.0 * trace(GammaA*GammaA*traceColour(q1*diquark));
-    //c = trace(q1*traceSpin(diquark)); //NO TRACESPIN???
+  //result = trace(q1*diquark); // TODO: Apologies, Felix - compiler errors
+  assert( 0 && "TODO: Felix, please fix prior line - compiler errors" );

    sliceSum(c,buf,Tp);

--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -25,7 +25,10 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Global.hpp>
+
+#include <chrono>
+#include <ctime>
+#include <Hadrons/Utilities/Contractor.hpp>
 #include <Hadrons/A2AMatrix.hpp>
 #include <Hadrons/DiskVector.hpp>
 #include <Hadrons/TimerArray.hpp>
@@ -34,66 +37,8 @@ using namespace Grid;
 using namespace QCD;
 using namespace Hadrons;

-#define TIME_MOD(t) (((t) + par.global.nt) % par.global.nt)
-
-namespace Contractor
-{
-    class TrajRange: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
-                                        unsigned int, start,
-                                        unsigned int, end,
-                                        unsigned int, step);
-    };
-    
-    class GlobalPar: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
-                                        TrajRange, trajCounter,
-                                        unsigned int, nt,
-                                        std::string, diskVectorDir,
-                                        std::string, output);
-    };
-
-    class A2AMatrixPar: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMatrixPar,
-                                        std::string, file,
-                                        std::string, dataset,
-                                        unsigned int, cacheSize,
-                                        std::string, name);
-    };
-
-    class ProductPar: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(ProductPar,
-                                        std::string, terms,
-                                        std::vector<std::string>, times,
-                                        std::string, translations,
-                                        bool, translationAverage);
-    };
-
-    class CorrelatorResult: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(CorrelatorResult,
-                                        std::vector<Contractor::A2AMatrixPar>,  a2aMatrix,
-                                        ProductPar, contraction,
-                                        std::vector<unsigned int>, times,
-                                        std::vector<ComplexD>, correlator);
-    };
-}
-
-struct ContractorPar
-{
-    Contractor::GlobalPar                  global;
-    std::vector<Contractor::A2AMatrixPar>  a2aMatrix;
-    std::vector<Contractor::ProductPar>    product;
-};
+// Separator to be used between contraction terms only (underscores elsewhere)
+std::string Separator{ "_" };

 void makeTimeSeq(std::vector<std::vector<unsigned int>> &timeSeq, 
                 const std::vector<std::set<unsigned int>> &times,
@@ -119,7 +64,7 @@ void makeTimeSeq(std::vector<std::vector<unsigned int>> &timeSeq,
 {
    std::vector<unsigned int> current(times.size());

-    makeTimeSeq(timeSeq, times, current, times.size());
+  makeTimeSeq(timeSeq, times, current, static_cast<unsigned int>(times.size()));
 }

 void saveCorrelator(const Contractor::CorrelatorResult &result, const std::string dir,
@@ -130,12 +75,12 @@ void saveCorrelator(const Contractor::CorrelatorResult &result, const std::strin

    for (unsigned int i = 0; i < terms.size() - 1; i++)
    {
-        fileStem += terms[i] + "_" + std::to_string(result.times[i]) + "_";
+        fileStem += terms[i] + Separator + std::to_string(result.times[i]) + Separator;
    }
    fileStem += terms.back();
    if (!result.contraction.translationAverage)
    {
-        fileStem += "_dt_" + std::to_string(dt);
+        fileStem += Separator + "dt" + Separator + std::to_string(dt);
    }
    filename = dir + "/" + RESULT_FILE_NAME(fileStem, traj);
    std::cout << "Saving correlator to '" << filename << "'" << std::endl;
@@ -240,30 +185,75 @@ int main(int argc, char* argv[])
 {
    // parse command line
    std::string parFilename;
+    bool        bOnlyWriteUsedA2AMatrices{ false };
+    int         ArgCount{ 0 };
+    bool        bCmdLineError{ false };
+    for( int i = 1; i < argc; i++ ) {
+        if( argv[i][0] == '-' ) {
+            // Switches
+            bool bSwitchOK = false;
+            switch( argv[i][1] ) {
+                case 'a':
+                    if( argv[i][2] == 0 ) {
+                        bOnlyWriteUsedA2AMatrices = true;
+                        bSwitchOK = true;
+                        std::cout << "Only A2AMatrices used in each contraction will be written" << std::endl;
+                    }
+                    break;
+                case 's':
+                    if( argv[i][2] )
+                        Separator = &argv[i][2];
+                    else
+                        Separator = ".";
+                    bSwitchOK = true;
+                    std::cout << "Using \"" << Separator << "\" as name separator" << std::endl;
+                    break;
+            }
+            if( !bSwitchOK ) {
+                std::cerr << "Urecognised switch \"" << argv[i] << "\"" << std::endl;
+                bCmdLineError = true;
+            }
+        } else {
+            // Arguments
+            switch( ++ArgCount ) {
+                case 1:
+                    parFilename = argv[i];
+                    break;
+                default:
+                    std::cerr << "Unused argument \"" << argv[i] << "\"" << std::endl;
+                    break;
+            }
+        }
+    }

-    if (argc != 2)
+    if (ArgCount != 1 or bCmdLineError)
    {
-        std::cerr << "usage: " << argv[0] << " <parameter file>";
-        std::cerr << std::endl;
+        std::cerr << "usage: " << argv[0] << " <parameter file>"
+                     "\n\t-a\tSimple Correlators (only describe A2AMatrices used for contraction)"
+                     "\n\t-s[sep]\tSeparator \"sep\" used between name components."
+                     "\n\t\tDefaults to \"_\", or \".\" if -s specified without sep"
+                  << std::endl;
        
        return EXIT_FAILURE;
    }
-    parFilename = argv[1];
+
+    // Log what file we're processing and when we started
+    const std::chrono::system_clock::time_point start{ std::chrono::system_clock::now() };
+    std::time_t now = std::chrono::system_clock::to_time_t( start );
+    std::cout << "Start " << parFilename << " " << std::ctime( &now );

    // parse parameter file
-    ContractorPar par;
-    unsigned int  nMat, nCont;
+    Contractor::ContractorPar par;
    XmlReader     reader(parFilename);

    read(reader, "global",    par.global);
    read(reader, "a2aMatrix", par.a2aMatrix);
    read(reader, "product",   par.product);
-    nMat  = par.a2aMatrix.size();
-    nCont = par.product.size();
+    const unsigned int nMat  { static_cast<unsigned int>(par.a2aMatrix.size()) };
+    const unsigned int nCont { static_cast<unsigned int>(par.product.size()) };

    // create diskvectors
    std::map<std::string, EigenDiskVector<ComplexD>> a2aMat;
-    unsigned int                                     cacheSize;

    for (auto &p: par.a2aMatrix)
    {
@@ -279,13 +269,15 @@ int main(int argc, char* argv[])
        std::cout << ":::::::: Trajectory " << traj << std::endl;

        // load data
+        int iSeq = 0;
        for (auto &p: par.a2aMatrix)
        {
            std::string filename = p.file;
-            double      t, size;
+            double      t;

            tokenReplace(filename, "traj", traj);
-            std::cout << "======== Loading '" << filename << "'" << std::endl;
+            std::cout << "======== Loading '" << filename << "'"
+                      << "\nA2AMatrix " << ++iSeq << " of " << nMat << " = " << p.name << std::endl;

            A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(filename, p.dataset, par.global.nt);

@@ -297,6 +289,7 @@ int main(int argc, char* argv[])
        // contract
        EigenDiskVector<ComplexD>::Matrix buf;

+        iSeq = 0;
        for (auto &p: par.product)
        {
            std::vector<std::string>               term = strToVec<std::string>(p.terms);
@@ -306,11 +299,11 @@ int main(int argc, char* argv[])
            std::vector<A2AMatrixTr<ComplexD>>     lastTerm(par.global.nt);
            A2AMatrix<ComplexD>                    prod, buf, tmp;
            TimerArray                             tAr;
-            double                                 fusec, busec, flops, bytes, tusec;
+            double                                 fusec, busec, flops, bytes;
            Contractor::CorrelatorResult           result;             

            tAr.startTimer("Total");
-            std::cout << "======== Contraction tr(";
+            std::cout << "======== Contraction " << ++iSeq << " of " << nCont << " tr(";
            for (unsigned int g = 0; g < term.size(); ++g)
            {
                std::cout << term[g] << ((g == term.size() - 1) ? ')' : '*');
@@ -328,7 +321,9 @@ int main(int argc, char* argv[])
            }
            for (auto &m: par.a2aMatrix)
            {
-                if (std::find(result.a2aMatrix.begin(), result.a2aMatrix.end(), m) == result.a2aMatrix.end())
+                // For simple correlators, only include A2AMatrix info for correlators in this contraction
+                if ( ( !bOnlyWriteUsedA2AMatrices or std::find( term.begin(), term.end(), m.name ) != term.end() )
+                  and std::find(result.a2aMatrix.begin(), result.a2aMatrix.end(), m) == result.a2aMatrix.end())
                {
                    result.a2aMatrix.push_back(m);
                    tokenReplace(result.a2aMatrix.back().file, "traj", traj);
@@ -451,5 +446,12 @@ int main(int argc, char* argv[])
        }
    }

+    // Mention that we're finished, what the time is and how long it took
+    const std::chrono::system_clock::time_point stop{ std::chrono::system_clock::now() };
+    now = std::chrono::system_clock::to_time_t( stop );
+    const std::chrono::duration<double> duration_seconds = stop - start;
+    const double hours{ ( duration_seconds.count() + 0.5 ) / 3600 };
+    std::cout << "Stop " << parFilename << " " << std::ctime( &now )
+              << "Total duration " << std::fixed << std::setprecision(1) << hours << " hours." << std::endl;
    return EXIT_SUCCESS;
 }
--- a/Hadrons/Utilities/Contractor.hpp
+++ b/Hadrons/Utilities/Contractor.hpp
@@ -36,4 +36,74 @@ BEGIN_HADRONS_NAMESPACE

 END_HADRONS_NAMESPACE

+#define BEGIN_CONTRACTOR_NAMESPACE namespace Contractor{
+BEGIN_CONTRACTOR_NAMESPACE
+
+using Grid::Serializable;
+using Grid::Reader;
+using Grid::Writer;
+using Grid::ComplexD;
+
+class TrajRange: Serializable
+{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
+                                  unsigned int, start,
+                                  unsigned int, end,
+                                  unsigned int, step);
+};
+
+class GlobalPar: Serializable
+{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
+                                  TrajRange, trajCounter,
+                                  unsigned int, nt,
+                                  std::string, diskVectorDir,
+                                  std::string, output);
+};
+
+class A2AMatrixPar: Serializable
+{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMatrixPar,
+                                  std::string, file,
+                                  std::string, dataset,
+                                  unsigned int, cacheSize,
+                                  std::string, name);
+};
+
+class ProductPar: Serializable
+{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(ProductPar,
+                                  std::string, terms,
+                                  std::vector<std::string>, times,
+                                  std::string, translations,
+                                  bool, translationAverage);
+};
+
+class CorrelatorResult: Serializable
+{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(CorrelatorResult,
+                                  std::vector<Contractor::A2AMatrixPar>,  a2aMatrix,
+                                  ProductPar, contraction,
+                                  std::vector<unsigned int>, times,
+                                  std::vector<ComplexD>, correlator);
+};
+
+struct ContractorPar
+{
+  Contractor::GlobalPar                  global;
+  std::vector<Contractor::A2AMatrixPar>  a2aMatrix;
+  std::vector<Contractor::ProductPar>    product;
+};
+
+// Useful ... so long as there's a ContractorPar named par in scope
+#define TIME_MOD(t) (((t) + par.global.nt) % par.global.nt)
+
+#define END_CONTRACTOR_NAMESPACE }
+END_CONTRACTOR_NAMESPACE
+
 #endif // Hadrons_Contractor_hpp_
--- a/Hadrons/modules.inc
+++ b/Hadrons/modules.inc
@@ -143,7 +143,6 @@ modules_hpp =\
  Modules/MContraction/WeakEye3pt.hpp \
  Modules/MContraction/WeakNonEye3pt.hpp \
  Modules/MContraction/Baryon.hpp \
-  Modules/MContraction/Baryon_old.hpp \
  Modules/MContraction/Meson.hpp \
  Modules/MContraction/A2ALoop.hpp \
  Modules/MContraction/Gamma3pt.hpp \