Hadrons: eigenpack support for multiprecision I/O

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 <typename F>
+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 <typename T, typename TIo = T>
+    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 <typename T, typename TIo = T>
+    static void readPack(std::vector<T> &evec, std::vector<RealD> &eval,
+                         PackRecord &record, const std::string filename, 
+                         const unsigned int size, bool multiFile, 
+                         GridBase *gridIo = nullptr)
+    {
+        std::unique_ptr<TIo> ioBuf{nullptr};
+        ScidacReader         binReader;
+
+        if (typeHash<T>() != typeHash<TIo>())
+        {
+            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 <typename T, typename TIo = T>
+    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 <typename T, typename TIo = T>
+    static void writePack(const std::string filename, std::vector<T> &evec, 
+                          std::vector<RealD> &eval, PackRecord &record, 
+                          const unsigned int size, bool multiFile, 
+                          GridBase *gridIo = nullptr)
+    {
+        GridBase             *grid = evec[0]._grid;
+        std::unique_ptr<TIo> ioBuf{nullptr}; 
+        std::unique_ptr<T>   testBuf{nullptr};
+        ScidacWriter         binWriter(grid->IsBoss());
+
+        if (typeHash<T>() != typeHash<TIo>())
+        {
+            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 <typename F, typename FIo = F>
 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<F>() != typeHash<FIo>())
+        {
+            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)
+        EigenPackIo::readPack<F, FIo>(evec, eval, record, evecFilename(fileStem, traj, multiFile), evec.size(), multiFile, gridIo_);
-        {
+        HADRONS_DUMP_EP_METADATA(record);
-            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;
-        }
     }
 
     virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
+        EigenPackIo::writePack<F, FIo>(evecFilename(fileStem, traj, multiFile), evec, eval, record, evec.size(), multiFile, gridIo_);
-        {
-            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 <typename T>
-    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 <typename T>
-    static void writeElement(ScidacWriter &binWriter, T &evec, VecRecord &vecRecord)
-    {
-        binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC);
     }
 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 <typename T>
+    
-    void basicRead(std::vector<T> &evec, std::vector<RealD> &eval,
+protected:
-                   const std::string filename, const unsigned int size)
+    GridBase *gridIo_;
-    {
-        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 <typename T>
-    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 <typename T>
-    void basicWrite(const std::string filename, std::vector<T> &evec, 
-                    const std::vector<RealD> &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 <typename T>
-    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();
-    }
 };
 
-template <typename FineF, typename CoarseF>
+template <typename FineF, typename CoarseF, 
-class CoarseEigenPack: public EigenPack<FineF>
+          typename FineFIo = FineF, typename CoarseFIo = CoarseF>
+class CoarseEigenPack: public EigenPack<FineF, FineFIo>
 {
 public:
-    typedef CoarseF CoarseField;
+    typedef CoarseF CoarseField;         
-public:
-    std::vector<RealD>   evalCoarse;
     std::vector<CoarseF> evecCoarse;
+    std::vector<RealD>   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<FineF>() != typeHash<FineFIo>())
+        {
+            if (gridFineIo == nullptr)
+            {
+                HADRONS_ERROR(Definition, 
+                              "Fine I/O type different from vector type but null fine I/O grid passed");
+            }
+        }
+        if (typeHash<CoarseF>() != typeHash<CoarseFIo>())
+        {
+            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<FineF>::resize(sizeFine, gridFine);
+        EigenPack<FineF, FineFIo>::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)
+        EigenPack<FineF, FineFIo>::read(fileStem + "_fine", multiFile, traj);
-        {
-            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());
-        }
     }
 
     void readCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
+        PackRecord dummy;
-        {
+
-            for(int k = 0; k < evecCoarse.size(); ++k)
+        EigenPackIo::readPack<CoarseF, CoarseFIo>(evecCoarse, evalCoarse, dummy, 
-            {
+                              this->evecFilename(fileStem + "_coarse", traj, multiFile), 
-                this->basicReadSingle(evecCoarse[k], evalCoarse[k], this->evecFilename(fileStem + "_coarse", k, traj), k);
+                              evecCoarse.size(), multiFile, gridCoarseIo_);
-            }
-        }
-        else
-        {
-            this->basicRead(evecCoarse, evalCoarse, this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse.size());
-        }
     }
 
     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)
+        EigenPack<FineF, FineFIo>::write(fileStem + "_fine", multiFile, traj);
-        {
-            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());
-        }
     }
 
     void writeCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
+        EigenPackIo::writePack<CoarseF, CoarseFIo>(this->evecFilename(fileStem + "_coarse", traj, multiFile), 
-        {
+                                                   evecCoarse, evalCoarse, this->record, 
-            for(int k = 0; k < evecCoarse.size(); ++k)
+                                                   evecCoarse.size(), multiFile, gridCoarseIo_);
-            {
-                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());
-        }
     }
     
     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 <typename FImpl>
+template <typename FImpl, typename FImplIo = FImpl>
-using FermionEigenPack = EigenPack<typename FImpl::FermionField>;
+using FermionEigenPack = EigenPack<typename FImpl::FermionField, typename FImplIo::FermionField>;
 
-template <typename FImpl, int nBasis>
+template <typename FImpl, int nBasis, typename FImplIo = FImpl>
 using CoarseFermionEigenPack = CoarseEigenPack<
     typename FImpl::FermionField,
     typename LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
                                    typename FImpl::SiteComplex, 
+                                   nBasis>::CoarseField,
+    typename FImplIo::FermionField,
+    typename LocalCoherenceLanczos<typename FImplIo::SiteSpinor, 
+                                   typename FImplIo::SiteComplex, 
                                    nBasis>::CoarseField>;
 
 #undef HADRONS_DUMP_EP_METADATA

Hadrons/Utilities/EigenPackCast.cc

@@ -59,8 +59,12 @@ void convert(const std::string outFilename, const std::string inFilename,
         }
     }
 
-    FOut bufOut(gOut);
+    FOut         bufOut(gOut);
-    FIn  bufIn(gIn), testIn(gIn);
+    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);
+            EigenPackIo::readHeader(record, binReader);
-            EPOut::writeHeader(binWriter, record);
+            EigenPackIo::writeHeader(binWriter, record);
-            EPIn::readElement(bufIn, vecRecord, binReader);
+            EigenPackIo::readElement<FIn>(bufIn, eval, k, binReader);
-            precisionChange(bufOut, bufIn);
+            EigenPackIo::writeElement<FIn, FOut>(binWriter, bufIn, eval, k, &bufOut, &testIn);
-            precisionChange(testIn, bufOut);
-            testIn -= bufIn;
-            LOG(Message) << "Diff norm^2: " << norm2(testIn) << std::endl;
-            EPOut::writeElement(binWriter, bufOut, vecRecord);
             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);
+        EigenPackIo::readHeader(record, binReader);
-        EPOut::writeHeader(binWriter, record);
+        EigenPackIo::writeHeader(binWriter, record);
         for(unsigned int k = 0; k < size; ++k)
         {
-            VecRecord vecRecord;
+            EigenPackIo::readElement<FIn>(bufIn, eval, k, binReader);
-
+            EigenPackIo::writeElement<FIn, FOut>(binWriter, bufIn, eval, k, &bufOut, &testIn);
-            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);
         }
         binWriter.close();
         binReader.close();