diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index b2a6ab53..b87d714b 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -143,8 +143,9 @@ void GlobalSharedMemory::OptimalCommunicator(const std::vector &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 &processors,Grid_MPI_Comm & optimal_comm) { diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 19fb03c5..b48a72c9 100644 --- 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::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 -//typename FImpl::PropagatorField BaryonUtils::quarkContract13(const PropagatorField &q1, -// const PropagatorField &q2) -void BaryonUtils::quarkContract13(const PropagatorField &q1, - const PropagatorField &q2, - PropagatorField &q_out) +LatticeSpinColourMatrix BaryonUtils::quarkContract13(const PropagatorField &q1, + const PropagatorField &q2) { GridBase *grid = q1._grid; - std::vector> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; std::vector 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;ssoSites();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::quarkContract13(const PropagatorField &q1, }}} } } + } //end loop over lattice sites - q_out._odata[ss]=D_out; - - } //end loop over lattice sites - - //return q_out; + return q_out; } - }} diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 047f0aff..85f265cf 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -60,7 +60,6 @@ #include #include #include -#include #include #include #include diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index f78759d8..7d44561e 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -116,6 +116,12 @@ template void TBaryon::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,26 +178,76 @@ void TBaryon::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); result.corr.resize(nt); const std::string gamma{ par().gamma }; std::vector 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 buf1; + std::vector buf2; + std::vector buf3; + std::vector buf4; + std::vector buf5; + std::vector buf6; const Gamma GammaA{ Gamma::Algebra::Identity }; const Gamma GammaB{ al }; - BaryonUtils::ContractBaryons(q1,q2,q3,GammaA,GammaB,c); + //BaryonUtils::ContractBaryons(q1,q2,q3,GammaA,GammaB,c); + BaryonUtils::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 diff --git a/Hadrons/Modules/MContraction/Baryon2.hpp b/Hadrons/Modules/MContraction/Baryon2.hpp index 7bf977cf..1d434d0a 100644 --- a/Hadrons/Modules/MContraction/Baryon2.hpp +++ b/Hadrons/Modules/MContraction/Baryon2.hpp @@ -182,11 +182,13 @@ void TBaryon2::execute(void) const Gamma GammaA{ Gamma::Algebra::Identity }; const Gamma GammaB{ al }; - //diquark = BaryonUtils::quarkContract13(q2*GammaB,GammaB*q3); - BaryonUtils::quarkContract13(q2*GammaB,GammaB*q3,diquark); + LatticeSpinColourMatrix diquark( q1._grid ); // TODO: Felix, I added "q1._grid". I presume this is correct? + + diquark = BaryonUtils::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); diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc index cfa4492e..7d643778 100644 --- 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 + +#include +#include +#include #include #include #include @@ -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, times, - std::string, translations, - bool, translationAverage); - }; - - class CorrelatorResult: Serializable - { - public: - GRID_SERIALIZABLE_CLASS_MEMBERS(CorrelatorResult, - std::vector, a2aMatrix, - ProductPar, contraction, - std::vector, times, - std::vector, correlator); - }; -} - -struct ContractorPar -{ - Contractor::GlobalPar global; - std::vector a2aMatrix; - std::vector product; -}; +// Separator to be used between contraction terms only (underscores elsewhere) +std::string Separator{ "_" }; void makeTimeSeq(std::vector> &timeSeq, const std::vector> ×, @@ -119,10 +64,10 @@ void makeTimeSeq(std::vector> &timeSeq, { std::vector current(times.size()); - makeTimeSeq(timeSeq, times, current, times.size()); + makeTimeSeq(timeSeq, times, current, static_cast(times.size())); } -void saveCorrelator(const Contractor::CorrelatorResult &result, const std::string dir, +void saveCorrelator(const Contractor::CorrelatorResult &result, const std::string dir, const unsigned int dt, const unsigned int traj) { std::string fileStem = "", filename; @@ -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; @@ -239,31 +184,76 @@ inline std::ostream & operator<< (std::ostream& s, const Bytes &&b) int main(int argc, char* argv[]) { // parse command line - std::string parFilename; + 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] << " "; - std::cerr << std::endl; + std::cerr << "usage: " << argv[0] << " " + "\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(par.a2aMatrix.size()) }; + const unsigned int nCont { static_cast(par.product.size()) }; // create diskvectors std::map> 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 a2aIo(filename, p.dataset, par.global.nt); @@ -297,6 +289,7 @@ int main(int argc, char* argv[]) // contract EigenDiskVector::Matrix buf; + iSeq = 0; for (auto &p: par.product) { std::vector term = strToVec(p.terms); @@ -306,11 +299,11 @@ int main(int argc, char* argv[]) std::vector> lastTerm(par.global.nt); A2AMatrix 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); @@ -450,6 +445,13 @@ int main(int argc, char* argv[]) printTimeProfile(tAr.getTimings(), tAr.getTimer("Total")); } } - + + // 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 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; } diff --git a/Hadrons/Utilities/Contractor.hpp b/Hadrons/Utilities/Contractor.hpp index decd13aa..27e55b50 100644 --- 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, times, + std::string, translations, + bool, translationAverage); +}; + +class CorrelatorResult: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(CorrelatorResult, + std::vector, a2aMatrix, + ProductPar, contraction, + std::vector, times, + std::vector, correlator); +}; + +struct ContractorPar +{ + Contractor::GlobalPar global; + std::vector a2aMatrix; + std::vector 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_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 36e6429b..5fd2f6d8 100644 --- 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 \