diff --git a/Hadrons/EigenPack.hpp b/Hadrons/EigenPack.hpp index 0be998b7..0e6f95d7 100644 --- a/Hadrons/EigenPack.hpp +++ b/Hadrons/EigenPack.hpp @@ -39,12 +39,12 @@ BEGIN_HADRONS_NAMESPACE #define HADRONS_DEFAULT_LANCZOS_NBASIS 60 #endif -#define HADRONS_DUMP_EP_METADATA \ +#define HADRONS_DUMP_EP_METADATA(record) \ LOG(Message) << "Eigenpack metadata:" << std::endl;\ LOG(Message) << "* operator" << std::endl;\ -LOG(Message) << record.operatorXml << std::endl;\ +LOG(Message) << (record).operatorXml << std::endl;\ LOG(Message) << "* solver" << std::endl;\ -LOG(Message) << record.solverXml << std::endl; +LOG(Message) << (record).solverXml << std::endl; struct PackRecord { @@ -59,7 +59,172 @@ struct VecRecord: Serializable VecRecord(void): index(0), eval(0.) {} }; -template +namespace EigenPackIo +{ + void readHeader(PackRecord &record, ScidacReader &binReader) + { + std::string recordXml; + + binReader.readLimeObject(recordXml, SCIDAC_FILE_XML); + XmlReader xmlReader(recordXml, true, "eigenPackPar"); + xmlReader.push(); + xmlReader.readCurrentSubtree(record.operatorXml); + xmlReader.nextElement(); + xmlReader.readCurrentSubtree(record.solverXml); + } + + template + void readElement(T &evec, RealD &eval, const unsigned int index, + ScidacReader &binReader, TIo *ioBuf = nullptr) + { + VecRecord vecRecord; + + LOG(Message) << "Reading eigenvector " << index << std::endl; + if (ioBuf == nullptr) + { + binReader.readScidacFieldRecord(evec, vecRecord); + } + else + { + binReader.readScidacFieldRecord(*ioBuf, vecRecord); + precisionChange(evec, *ioBuf); + } + if (vecRecord.index != index) + { + HADRONS_ERROR(Io, "Eigenvector " + std::to_string(index) + " has a" + + " wrong index (expected " + std::to_string(vecRecord.index) + + ")"); + } + eval = vecRecord.eval; + } + + template + static void readPack(std::vector &evec, std::vector &eval, + PackRecord &record, const std::string filename, + const unsigned int size, bool multiFile, + GridBase *gridIo = nullptr) + { + std::unique_ptr ioBuf{nullptr}; + ScidacReader binReader; + + if (typeHash() != typeHash()) + { + if (gridIo == nullptr) + { + HADRONS_ERROR(Definition, + "I/O type different from vector type but null I/O grid passed"); + } + ioBuf.reset(new TIo(gridIo)); + } + if (multiFile) + { + std::string fullFilename; + + for(int k = 0; k < size; ++k) + { + fullFilename = filename + "/v" + std::to_string(k) + ".bin"; + binReader.open(fullFilename); + readHeader(record, binReader); + readElement(evec[k], eval[k], k, binReader, ioBuf.get()); + binReader.close(); + } + } + else + { + binReader.open(filename); + readHeader(record, binReader); + for(int k = 0; k < size; ++k) + { + readElement(evec[k], eval[k], k, binReader, ioBuf.get()); + } + binReader.close(); + } + } + + void writeHeader(ScidacWriter &binWriter, PackRecord &record) + { + XmlWriter xmlWriter("", "eigenPackPar"); + + xmlWriter.pushXmlString(record.operatorXml); + xmlWriter.pushXmlString(record.solverXml); + binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML); + } + + template + void writeElement(ScidacWriter &binWriter, T &evec, RealD &eval, + const unsigned int index, TIo *ioBuf, + T *testBuf = nullptr) + { + VecRecord vecRecord; + + LOG(Message) << "Writing eigenvector " << index << std::endl; + vecRecord.eval = eval; + vecRecord.index = index; + if ((ioBuf == nullptr) || (testBuf == nullptr)) + { + binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC); + } + else + { + precisionChange(*ioBuf, evec); + precisionChange(*testBuf, *ioBuf); + *testBuf -= evec; + LOG(Message) << "Precision diff norm^2 " << norm2(*testBuf) << std::endl; + binWriter.writeScidacFieldRecord(*ioBuf, vecRecord, DEFAULT_ASCII_PREC); + } + } + + template + static void writePack(const std::string filename, std::vector &evec, + std::vector &eval, PackRecord &record, + const unsigned int size, bool multiFile, + GridBase *gridIo = nullptr) + { + GridBase *grid = evec[0]._grid; + std::unique_ptr ioBuf{nullptr}; + std::unique_ptr testBuf{nullptr}; + ScidacWriter binWriter(grid->IsBoss()); + + if (typeHash() != typeHash()) + { + if (gridIo == nullptr) + { + HADRONS_ERROR(Definition, + "I/O type different from vector type but null I/O grid passed"); + } + ioBuf.reset(new TIo(gridIo)); + testBuf.reset(new T(grid)); + } + if (multiFile) + { + std::string fullFilename; + + for(int k = 0; k < size; ++k) + { + fullFilename = filename + "/v" + std::to_string(k) + ".bin"; + + makeFileDir(fullFilename, grid); + binWriter.open(fullFilename); + writeHeader(binWriter, record); + writeElement(binWriter, evec[k], eval[k], k, ioBuf.get(), testBuf.get()); + binWriter.close(); + } + } + else + { + makeFileDir(filename, grid); + binWriter.open(filename); + writeHeader(binWriter, record); + for(int k = 0; k < size; ++k) + { + writeElement(binWriter, evec[k], eval[k], k, ioBuf.get(), testBuf.get()); + } + binWriter.close(); + } + } +} + +template class EigenPack { public: @@ -72,8 +237,17 @@ public: EigenPack(void) = default; virtual ~EigenPack(void) = default; - EigenPack(const size_t size, GridBase *grid) + EigenPack(const size_t size, GridBase *grid, GridBase *gridIo = nullptr) { + if (typeHash() != typeHash()) + { + if (gridIo == nullptr) + { + HADRONS_ERROR(Definition, + "I/O type different from vector type but null I/O grid passed"); + } + } + gridIo_ = gridIo; resize(size, grid); } @@ -85,225 +259,92 @@ public: virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1) { - if (multiFile) - { - for(int k = 0; k < evec.size(); ++k) - { - basicReadSingle(evec[k], eval[k], evecFilename(fileStem, k, traj), k); - if (k == 0) - { - HADRONS_DUMP_EP_METADATA; - } - } - } - else - { - basicRead(evec, eval, evecFilename(fileStem, -1, traj), evec.size()); - HADRONS_DUMP_EP_METADATA; - } + EigenPackIo::readPack(evec, eval, record, evecFilename(fileStem, traj, multiFile), evec.size(), multiFile, gridIo_); + HADRONS_DUMP_EP_METADATA(record); } virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1) { - if (multiFile) - { - for(int k = 0; k < evec.size(); ++k) - { - basicWriteSingle(evecFilename(fileStem, k, traj), evec[k], eval[k], k); - } - } - else - { - basicWrite(evecFilename(fileStem, -1, traj), evec, eval, evec.size()); - } - } - - static void readHeader(PackRecord &record, ScidacReader &binReader) - { - std::string recordXml; - - binReader.readLimeObject(recordXml, SCIDAC_FILE_XML); - XmlReader xmlReader(recordXml, true, "eigenPackPar"); - xmlReader.push(); - xmlReader.readCurrentSubtree(record.operatorXml); - xmlReader.nextElement(); - xmlReader.readCurrentSubtree(record.solverXml); - } - - template - static void readElement(T &evec, VecRecord &vecRecord, ScidacReader &binReader) - { - binReader.readScidacFieldRecord(evec, vecRecord); - } - - static void writeHeader(ScidacWriter &binWriter, PackRecord &record) - { - XmlWriter xmlWriter("", "eigenPackPar"); - - xmlWriter.pushXmlString(record.operatorXml); - xmlWriter.pushXmlString(record.solverXml); - binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML); - } - - template - static void writeElement(ScidacWriter &binWriter, T &evec, VecRecord &vecRecord) - { - binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC); + EigenPackIo::writePack(evecFilename(fileStem, traj, multiFile), evec, eval, record, evec.size(), multiFile, gridIo_); } protected: - std::string evecFilename(const std::string stem, const int vec, const int traj) + std::string evecFilename(const std::string stem, const int traj, const bool multiFile) { std::string t = (traj < 0) ? "" : ("." + std::to_string(traj)); - if (vec == -1) + if (multiFile) { - return stem + t + ".bin"; + return stem + t; } else { - return stem + t + "/v" + std::to_string(vec) + ".bin"; + return stem + t + ".bin"; } } - template - void basicRead(std::vector &evec, std::vector &eval, - const std::string filename, const unsigned int size) - { - ScidacReader binReader; - - binReader.open(filename); - readHeader(record, binReader); - for(int k = 0; k < size; ++k) - { - VecRecord vecRecord; - - LOG(Message) << "Reading eigenvector " << k << std::endl; - readElement(evec[k], vecRecord, binReader); - if (vecRecord.index != k) - { - HADRONS_ERROR(Io, "Eigenvector " + std::to_string(k) + " has a" - + " wrong index (expected " + std::to_string(vecRecord.index) - + ") in file '" + filename + "'"); - } - eval[k] = vecRecord.eval; - } - binReader.close(); - } - - template - void basicReadSingle(T &evec, RealD &eval, const std::string filename, - const unsigned int index) - { - ScidacReader binReader; - VecRecord vecRecord; - - binReader.open(filename); - readHeader(record, binReader); - LOG(Message) << "Reading eigenvector " << index << std::endl; - readElement(evec, vecRecord, binReader); - if (vecRecord.index != index) - { - HADRONS_ERROR(Io, "Eigenvector " + std::to_string(index) + " has a" - + " wrong index (expected " + std::to_string(vecRecord.index) - + ") in file '" + filename + "'"); - } - eval = vecRecord.eval; - binReader.close(); - } - - template - void basicWrite(const std::string filename, std::vector &evec, - const std::vector &eval, const unsigned int size) - { - ScidacWriter binWriter(evec[0]._grid->IsBoss()); - - makeFileDir(filename, evec[0]._grid); - binWriter.open(filename); - writeHeader(binWriter, record); - for(int k = 0; k < size; ++k) - { - VecRecord vecRecord; - - vecRecord.index = k; - vecRecord.eval = eval[k]; - LOG(Message) << "Writing eigenvector " << k << std::endl; - writeElement(binWriter, evec[k], vecRecord); - } - binWriter.close(); - } - - template - void basicWriteSingle(const std::string filename, T &evec, - const RealD eval, const unsigned int index) - { - ScidacWriter binWriter(evec._grid->IsBoss()); - VecRecord vecRecord; - - makeFileDir(filename, evec._grid); - binWriter.open(filename); - writeHeader(binWriter, record); - vecRecord.index = index; - vecRecord.eval = eval; - LOG(Message) << "Writing eigenvector " << index << std::endl; - writeElement(binWriter, evec, vecRecord); - binWriter.close(); - } + +protected: + GridBase *gridIo_; }; -template -class CoarseEigenPack: public EigenPack +template +class CoarseEigenPack: public EigenPack { public: - typedef CoarseF CoarseField; -public: - std::vector evalCoarse; + typedef CoarseF CoarseField; std::vector evecCoarse; + std::vector evalCoarse; public: CoarseEigenPack(void) = default; virtual ~CoarseEigenPack(void) = default; CoarseEigenPack(const size_t sizeFine, const size_t sizeCoarse, - GridBase *gridFine, GridBase *gridCoarse) + GridBase *gridFine, GridBase *gridCoarse, + GridBase *gridFineIo = nullptr, + GridBase *gridCoarseIo = nullptr) { + if (typeHash() != typeHash()) + { + if (gridFineIo == nullptr) + { + HADRONS_ERROR(Definition, + "Fine I/O type different from vector type but null fine I/O grid passed"); + } + } + if (typeHash() != typeHash()) + { + if (gridCoarseIo == nullptr) + { + HADRONS_ERROR(Definition, + "Coarse I/O type different from vector type but null coarse I/O grid passed"); + } + } + this->gridIo_ = gridFineIo; + gridCoarseIo_ = gridCoarseIo; resize(sizeFine, sizeCoarse, gridFine, gridCoarse); } void resize(const size_t sizeFine, const size_t sizeCoarse, GridBase *gridFine, GridBase *gridCoarse) { - EigenPack::resize(sizeFine, gridFine); + EigenPack::resize(sizeFine, gridFine); evalCoarse.resize(sizeCoarse); evecCoarse.resize(sizeCoarse, gridCoarse); } void readFine(const std::string fileStem, const bool multiFile, const int traj = -1) { - if (multiFile) - { - for(int k = 0; k < this->evec.size(); ++k) - { - this->basicReadSingle(this->evec[k], this->eval[k], this->evecFilename(fileStem + "_fine", k, traj), k); - } - } - else - { - this->basicRead(this->evec, this->eval, this->evecFilename(fileStem + "_fine", -1, traj), this->evec.size()); - } + EigenPack::read(fileStem + "_fine", multiFile, traj); } void readCoarse(const std::string fileStem, const bool multiFile, const int traj = -1) { - if (multiFile) - { - for(int k = 0; k < evecCoarse.size(); ++k) - { - this->basicReadSingle(evecCoarse[k], evalCoarse[k], this->evecFilename(fileStem + "_coarse", k, traj), k); - } - } - else - { - this->basicRead(evecCoarse, evalCoarse, this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse.size()); - } + PackRecord dummy; + + EigenPackIo::readPack(evecCoarse, evalCoarse, dummy, + this->evecFilename(fileStem + "_coarse", traj, multiFile), + evecCoarse.size(), multiFile, gridCoarseIo_); } virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1) @@ -314,32 +355,14 @@ public: void writeFine(const std::string fileStem, const bool multiFile, const int traj = -1) { - if (multiFile) - { - for(int k = 0; k < this->evec.size(); ++k) - { - this->basicWriteSingle(this->evecFilename(fileStem + "_fine", k, traj), this->evec[k], this->eval[k], k); - } - } - else - { - this->basicWrite(this->evecFilename(fileStem + "_fine", -1, traj), this->evec, this->eval, this->evec.size()); - } + EigenPack::write(fileStem + "_fine", multiFile, traj); } void writeCoarse(const std::string fileStem, const bool multiFile, const int traj = -1) { - if (multiFile) - { - for(int k = 0; k < evecCoarse.size(); ++k) - { - this->basicWriteSingle(this->evecFilename(fileStem + "_coarse", k, traj), evecCoarse[k], evalCoarse[k], k); - } - } - else - { - this->basicWrite(this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse, evalCoarse, evecCoarse.size()); - } + EigenPackIo::writePack(this->evecFilename(fileStem + "_coarse", traj, multiFile), + evecCoarse, evalCoarse, this->record, + evecCoarse.size(), multiFile, gridCoarseIo_); } virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1) @@ -347,16 +370,22 @@ public: writeFine(fileStem, multiFile, traj); writeCoarse(fileStem, multiFile, traj); } +private: + GridBase *gridCoarseIo_; }; -template -using FermionEigenPack = EigenPack; +template +using FermionEigenPack = EigenPack; -template +template using CoarseFermionEigenPack = CoarseEigenPack< typename FImpl::FermionField, typename LocalCoherenceLanczos::CoarseField, + typename FImplIo::FermionField, + typename LocalCoherenceLanczos::CoarseField>; #undef HADRONS_DUMP_EP_METADATA diff --git a/Hadrons/Utilities/EigenPackCast.cc b/Hadrons/Utilities/EigenPackCast.cc index cf35a550..7247c92d 100644 --- a/Hadrons/Utilities/EigenPackCast.cc +++ b/Hadrons/Utilities/EigenPackCast.cc @@ -59,8 +59,12 @@ void convert(const std::string outFilename, const std::string inFilename, } } - FOut bufOut(gOut); - FIn bufIn(gIn), testIn(gIn); + FOut bufOut(gOut); + FIn bufIn(gIn), testIn(gIn); + ScidacWriter binWriter(gOut->IsBoss()); + ScidacReader binReader; + PackRecord record; + RealD eval; LOG(Message) << "==== EIGENPACK CONVERSION" << std::endl; LOG(Message) << "Lattice : " << gIn->GlobalDimensions() << std::endl; @@ -75,10 +79,6 @@ void convert(const std::string outFilename, const std::string inFilename, { for(unsigned int k = 0; k < size; ++k) { - ScidacWriter binWriter(gOut->IsBoss()); - ScidacReader binReader; - PackRecord record; - VecRecord vecRecord; std::string outV = outFilename + "/v" + std::to_string(k) + ".bin"; std::string inV = inFilename + "/v" + std::to_string(k) + ".bin"; @@ -88,40 +88,25 @@ void convert(const std::string outFilename, const std::string inFilename, makeFileDir(outV, gOut); binWriter.open(outV); binReader.open(inV); - EPIn::readHeader(record, binReader); - EPOut::writeHeader(binWriter, record); - EPIn::readElement(bufIn, vecRecord, binReader); - precisionChange(bufOut, bufIn); - precisionChange(testIn, bufOut); - testIn -= bufIn; - LOG(Message) << "Diff norm^2: " << norm2(testIn) << std::endl; - EPOut::writeElement(binWriter, bufOut, vecRecord); + EigenPackIo::readHeader(record, binReader); + EigenPackIo::writeHeader(binWriter, record); + EigenPackIo::readElement(bufIn, eval, k, binReader); + EigenPackIo::writeElement(binWriter, bufIn, eval, k, &bufOut, &testIn); binWriter.close(); binReader.close(); } } else { - ScidacWriter binWriter(gOut->IsBoss()); - ScidacReader binReader; - PackRecord record; - makeFileDir(outFilename, gOut); binWriter.open(outFilename); binReader.open(inFilename); - EPIn::readHeader(record, binReader); - EPOut::writeHeader(binWriter, record); + EigenPackIo::readHeader(record, binReader); + EigenPackIo::writeHeader(binWriter, record); for(unsigned int k = 0; k < size; ++k) { - VecRecord vecRecord; - - LOG(Message) << "==== Converting vector " << k << std::endl; - EPIn::readElement(bufIn, vecRecord, binReader); - precisionChange(bufOut, bufIn); - precisionChange(testIn, bufOut); - testIn -= bufIn; - LOG(Message) << "Diff norm^2: " << norm2(testIn) << std::endl; - EPOut::writeElement(binWriter, bufOut, vecRecord); + EigenPackIo::readElement(bufIn, eval, k, binReader); + EigenPackIo::writeElement(binWriter, bufIn, eval, k, &bufOut, &testIn); } binWriter.close(); binReader.close();