From c69a3b6ef6cedb173e7a56908d48c5850146c263 Mon Sep 17 00:00:00 2001 From: Michael Marshall <43034299+mmphys@users.noreply.github.com> Date: Wed, 29 Jan 2020 21:20:20 +0000 Subject: [PATCH] When saving eigenvectors, LapEvec now saves eigenvalues for every timeslice as well. I.e. nT x nVec eigenvalues are saved in FileName.evals.conf.h5. A new named tensor, "TimesliceEvals" can be used to simplify restoring these from disk. NB: The changes in BaseIO add support so that Eigen tensors can be easily used in MPI operations, e.g. GlobalSum. See LapEvec.hpp for an example of how this is done. --- Grid/serialisation/BaseIO.h | 36 ++++++++++++------------- Hadrons/Modules/MDistil/LapEvec.hpp | 23 ++++++++++++++-- Hadrons/Modules/MDistil/Perambulator.cc | 3 +++ Hadrons/NamedTensor.hpp | 14 ++++++++++ 4 files changed, 56 insertions(+), 20 deletions(-) diff --git a/Grid/serialisation/BaseIO.h b/Grid/serialisation/BaseIO.h index b858c835..bf424fc7 100644 --- a/Grid/serialisation/BaseIO.h +++ b/Grid/serialisation/BaseIO.h @@ -97,6 +97,23 @@ namespace Grid { template struct is_tensor_variable : public std::false_type {}; template struct is_tensor_variable::value && !is_tensor_fixed::value>::type> : public std::true_type {}; + + // Helper functions to get the ultimate scalar inside a tensor, and corresponding size + template + inline typename std::enable_if::value, const typename ET::Index>::type + getScalarCount(const ET &eigenTensor) { return eigenTensor.size() * Traits::count; } + template + inline typename std::enable_if::value, const typename ET::Scalar *>::type + getFirstScalar(const ET &eigenTensor) { return eigenTensor.data(); } + template + inline typename std::enable_if::value, typename ET::Scalar *>::type + getFirstScalar(ET &eigenTensor) { return eigenTensor.data(); } + template + inline typename std::enable_if::value, const typename Traits::scalar_type *>::type + getFirstScalar(const ET &eigenTensor) { return eigenTensor.data()->begin(); } + template + inline typename std::enable_if::value, typename Traits::scalar_type *>::type + getFirstScalar(ET &eigenTensor) { return eigenTensor.data()->begin(); } } // Abstract writer/reader classes //////////////////////////////////////////// @@ -128,23 +145,6 @@ namespace Grid { typename std::enable_if::value>::type write(const std::string &s, const ETensor &output); - // Helper functions for Scalar vs Container specialisations - template - inline typename std::enable_if::value, - const typename ETensor::Scalar *>::type - getFirstScalar(const ETensor &output) - { - return output.data(); - } - - template - inline typename std::enable_if::value, - const typename EigenIO::Traits::scalar_type *>::type - getFirstScalar(const ETensor &output) - { - return output.data()->begin(); - } - template inline typename std::enable_if::value, void>::type copyScalars(S * &pCopy, const S &Source) @@ -323,7 +323,7 @@ namespace Grid { std::vector CopyBuffer; const Index TotalNumElements = NumElements * Traits::count; if( !CopyData ) { - pWriteBuffer = getFirstScalar( output ); + pWriteBuffer = EigenIO::getFirstScalar( output ); } else { // Regardless of the Eigen::Tensor storage order, the copy will be Row Major CopyBuffer.resize( TotalNumElements ); diff --git a/Hadrons/Modules/MDistil/LapEvec.hpp b/Hadrons/Modules/MDistil/LapEvec.hpp index 9b5197c2..3c1122ca 100644 --- a/Hadrons/Modules/MDistil/LapEvec.hpp +++ b/Hadrons/Modules/MDistil/LapEvec.hpp @@ -263,6 +263,11 @@ void TLapEvec::execute(void) const int Ntlocal{gridHD->LocalDimensions()[Tdir]}; const int Ntfirst{gridHD->LocalStarts()[Tdir]}; uint32_t ConvergenceErrors{0}; + const int NtFull{env().getDim(Tdir)}; + TimesliceEvals Evals{ NtFull, LPar.Nvec }; + for (int t = 0; t < NtFull; t++) + for (int v = 0; v < LPar.Nvec; v++) + Evals.tensor( t, v ) = 0; for (int t = 0; t < Ntlocal; t++ ) { LOG(Message) << "------------------------------------------------------------" << std::endl; @@ -302,6 +307,8 @@ void TLapEvec::execute(void) InsertSliceLocal(eig[t].evec[i],eig4d.evec[i],0,t,Tdir); if(t==0 && Ntfirst==0) eig4d.eval[i] = eig[t].eval[i]; // TODO: Discuss: is this needed? Is there a better way? + if(gridLD->IsBoss()) // Only do this on one node per timeslice, so a global sum will work + Evals.tensor(t + Ntfirst,i) = eig[t].eval[i]; } } GridLogIRL.Active( PreviousIRLLogState ); @@ -322,9 +329,21 @@ void TLapEvec::execute(void) sOperatorXml.append( b->parString() ); sOperatorXml.append( "" ); eig4d.record.operatorXml = sOperatorXml; - sEigenPackName.append("."); - sEigenPackName.append(std::to_string(vm().getTrajectory())); + sEigenPackName.append(1, '.'); + std::size_t NameLen{ sEigenPackName.length() }; + const std::string sTrajNum{std::to_string(vm().getTrajectory())}; + sEigenPackName.append(sTrajNum); eig4d.write(sEigenPackName,false); + // Communicate eig[t].evec to boss-node, save into new object evecs + gridHD->GlobalSumVector(EigenIO::getFirstScalar(Evals.tensor), + static_cast(EigenIO::getScalarCount(Evals.tensor))); + if(gridHD->IsBoss()) + { + sEigenPackName.resize(NameLen); + sEigenPackName.append("evals."); + sEigenPackName.append(sTrajNum); + Evals.write( sEigenPackName ); + } } } diff --git a/Hadrons/Modules/MDistil/Perambulator.cc b/Hadrons/Modules/MDistil/Perambulator.cc index e710698e..a1b24a2d 100644 --- a/Hadrons/Modules/MDistil/Perambulator.cc +++ b/Hadrons/Modules/MDistil/Perambulator.cc @@ -53,5 +53,8 @@ const std::array NoiseTensor::DefaultIndexNames__{"nNoise", "nT" const std::string PerambTensor::Name__{"Perambulator"}; const std::array PerambTensor::DefaultIndexNames__{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"}; +const std::string TimesliceEvals::Name__{"TimesliceEigenValues"}; +const std::array TimesliceEvals::DefaultIndexNames__{"nT", "nVec"}; + END_MODULE_NAMESPACE END_HADRONS_NAMESPACE diff --git a/Hadrons/NamedTensor.hpp b/Hadrons/NamedTensor.hpp index 003a92fc..53e10e04 100644 --- a/Hadrons/NamedTensor.hpp +++ b/Hadrons/NamedTensor.hpp @@ -185,6 +185,20 @@ class PerambTensor : public NamedTensor : NamedTensor{Name__, DefaultIndexNames__, nT, nVec, LI, nNoise, nT_inv, SI} {} }; +class TimesliceEvals : public NamedTensor +{ + public: + static const std::string Name__; + static const std::array DefaultIndexNames__; + // Default constructor (assumes tensor will be loaded from file) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TimesliceEvals() : NamedTensor{Name__, DefaultIndexNames__} {} + + // Construct a named tensor explicitly specifying size of each dimension + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TimesliceEvals(Eigen::Index nT, Eigen::Index nVec) + : NamedTensor{Name__, DefaultIndexNames__, nT, nVec} {} +}; + END_MODULE_NAMESPACE END_HADRONS_NAMESPACE #endif // Hadrons_NamedTensor_hpp_